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ABSTRACT 

The  equations  of  motion  of  a  flexible  shut 1 1  e-beam- ant enna 
system  are  developed  and  discretized  using  an  assumed  modes 
approximation.  The  system  was  modeled  as  a  cantilever  beam 
rigidly  attached  to  the  shuttle  with  a  rigid  antenna  attached 
to  the  free  end  of  the  beam.  The  mass  and  dimension  data  for 
the  model  was  taken  from  a  NASA/IEEE  Design  Challenge  Paper 
[2]  dated  June  1984.  The  equations  of  motion  for  both  the 
shuttle-beam-antenna  rigid  body  movement  and  the  vibration  of 
the  beam  with  respect  to  the  shuttle  were  developed  making 
some  simplifyng  modifications  to  fit  the  modeling  as¬ 
sumptions.  Two  proof-mass  actuators,  capable  of  producing  a 
force  in  the  x  and  y  directions  only  with  no  torsional  con¬ 
trol  about  the  z  axis,  were  modeled  at  positions  along  the 
beam.  The  moments  resulting  from  any  torque  on  the  shuttle 
(due  to  reaction  jets  firing,  for  example),  and  moments 
applied  to  the  antenna  at  the  attach  point  were  also 
modeled.  The  equations  of  motion,  with  the  forces  and 
moments  evaluated,  were  put  in  matrix  form.  The  matrices 
were  diagonalized,  resulting  in  an  identity  mass  matrix  and 
diagonal  damping  and  stiffness  matrices. 

A  controller  was  developed  for  a  cursory  investigation 
into  the  controllability  of  the  system.  The  development  made 
use  of  linear  optimal  regulator  techniques  which  produce 
feedback  gains  proportional  to  the  state.  The  state  was 


truncated  to  the  amplitudes  and  velocities  of  the  twelve 
modes  having  the  lowest  frequencies.  Since  these  amplitudes 
and  velocities  could  not  be  measured  directly,  state  estima¬ 
tion  was  used.  The  feedback  gains  were  developed  using 
steady  state  optimal  regulator  theory.  The  closed  loop 
damping  coefficient  was  used  as  a  measure  of  control  improve¬ 
ment.  The  system  was  shown  to  be  stable  on  the  very  first 
control  attempt,  with  a  closed  loop  damping  coefficient  bet¬ 
ter  than  the  targetted  value.  Elimination  of  observation 
spillover  improved  the  controllability  slightly  for  this 
first  run.  More  runs  were  made  with  different  weightings  of 
the  controlled  modes  with  similar  results. 
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LARGE  SPACE  STRUCTURE  AS  APPLIED  TO 
A  SHUTTLE- ANTENNA  CONFIGURATION 

I .  Introduct ion 

With  the  advent  of  the  Space  Shuttle,  the  opportunities 
to  place  large,  flexible  strnctnres  into  space  are  becoming 
more  and  more  commonplace.  These  large  strnctnres,  along 
with  all  their  advantages,  bring  with  them  control  problems 
on  a  scale  never  before  encountered.  Extensive  research  is 
being  conducted  to  solve  these  control  problems — research 
which  is  constantly  being  updated  and  improved.  The  problem 
of  controlling  a  large  structure  containing  a  virtually  in 
finite  number  of  vibrational  modes  with  limited  onboard  com¬ 
puter  resources,  sensors,  and  actuators  has  become  the  focus 
for  intense  study.  Of  the  control  techniques  attempted  to 
date,  modern  state  space  control  methods  seem  to  be  the  most 
promising. 

In  January  1984  an  NASA  Design  Challenge  [1]  was  of¬ 
fered.  It  is  called  the  Spacecraft  Control  Laboratory  Ex¬ 
periment  (SCOLE)  and  serves  as  the  focus  of  a  design  chal¬ 
lenge  for  the  purpose  of  comparing  different  approaches  to 
control  synthesis,  modeling,  order  reduction,  state  estima¬ 
tion,  and  system  identification.  The  SCOLE  itself  is  pre- 
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sented  as  a  large  antenna  attached  to  the  shuttle  by  a  flex¬ 
ible  beam.  The  Challenge  consisted  of  two  parts — -  a  mathema¬ 
tical  analysis  and  a  laboratory  experiment.  Only  the  mathe¬ 
matical  analysis  is  addressed  in  this  paper. 

The  first  portion  of  this  paper  concerns  the  mathematical 
modeling  of  the  shut tl e-b e am- ante nna  system.  Unfortunately, 
there  were  many  errors  in  the  reference  [1].  The  paper  was 
re-released  in  June  1984  [2]  with  some  of  the  errors  cor¬ 
rected.  There  were,  however,  still  numerous  errors  in  [2] 
requiring  that  almost  all  of  the  mathematical  modeling  be 
done  from  scratch,  precluding  an  extensive  investigation  into 
the  controlling  of  the  system.  This  thesis  therefore  pro¬ 
vides  a  detailed  mathematical  model  of  the  system.  The 
system  equations  of  motion  were  developed  assuming  the  shut¬ 
tle  and  antenna  to  be  rigid  and  the  beam  to  be  flexible.  The 
beam  was  assumed  to  be  capable  of  transverse  bending  in  each 
of  two  orthogonal  directions  and  to  undergo  torsional  motion 
about  it's  long  axis.  The  assumed  modes  method  was  used  to 
discretize  the  beam  motion  and  a  set  of  linear  first-order 
equations  were  developed  for  the  system. 

Active  control  of  the  system  was  achieved  using  a  trun¬ 
cated  dynamical  model,  using  linear  optimal  regulator  theory 
and  modal  suppression  techniques  as  outlined  in  [6]  and  [7], 


VJil 


The  physical  model  of  the  SCOLE  is  shown  in  Fig  1.  It 
consists  of  the  shuttle,  a  130  foot  flexible  beam  attached  to 
the  shuttle's  center  of  mass  (an  assumption  made  for  modeling 
purposes),  and  a  rigid  antenna  attached  at  one  corner  to  the 
beam.  The  axis  system  is  as  shown:  The  x  (or  roll)  axis 
points  out  of  the  nose  of  the  shuttle,  the  y  (or  pitch)  axis 
points  out  the  shuttle's  right  wing,  and  the  z  (or  yaw)  axis 
points  out  the  bottom  of  the  shuttle,  which  is  nominally 
toward  the  Earth.  The  xyz  reference  frame  is  considered 
attached  to  the  shuttle  with  its  origin  at  the  center  of  mass 
(beam  attachment  point).  This  frame  is  f r?e  to  rotate  about 
an  XTZ  frame,  in  which  the  X  axis  points  along  the  velocity 
vector  of  the  shuttle  in  orbit,  the  Z  axis  points  to  the 
Earth's  center,  and  the  T  axis  completes  the  right-handed 
system.  The  XYZ  frame  is  considered  inertial  for  the  pur¬ 
poses  of  this  paper. 

Another  axis  system,  fixed  to  the  antenna  at  the  beam- 
antenna  attachment  point,  is  shown  in  Fig  2.  When  the  beam 
is  not  deformed,  this  x^y^z^  axis  system  aligns  itself  with 
the  xyz  axis  system.  When  the  beam  is  deformed,  the  *474*4 
axes  will  be  displaced  from  xyz  by  what  are  assumed  to  be 
very  small  angles.  This  axis  system  will  be  used  in  the 


mathematical  development  of  the  equations  of  motion  in  the 
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next  chapter. 

Reference  [2]  furnished  the  information  on  the  values  for 
the  masses  and  moments  of  inertia  for  the  shuttle,  reflector, 
and  entire  shnt tl e-be am- ant enna  system.  This  information  is 
presented  in  Table  I.  Also  included  are  the  modal  damping 
and  stiffness  coefficients  for  the  vibration  of  the  beam, 
along  with  the  masses  of  the  proof-mass  actuators. 

Before  developing  the  mathematical  model,  the  main  types 
of  motion  will  be  described,  along  with  the  assumptions  made 
for  both  simplification  and  clarification. 

Referring  to  Fig  1,  the  first  type  of  motion  considered 
is  the  tumbling  of  the  entire  sh u 1 1 1 e-b e am - an t e nna  system, 
expressed  as  some  arbitrary  rotation  of  the  xyz  frame  about 
the  inertial  XYZ  frame.  The  angular  velocity  is  given  as  u 
(no  subscript).  A  fundamental  assumption  for  the  purposes  of 
this  paper  is  that  the  shuttle  wishes  to  stay  aligned  with 
the  XTZ  frame,  so  any  rotation  out  of  that  alignment  will  be 
met  with  the  shuttle  firing  its  attitude  control  jets.  This 
means  that  any  rotation  will  be  small  and  of  short  duration. 
The  firing  of  the  reaction  jets  will  cause  the  beam  to  flex, 
resulting  in  the  second  type  of  motion  to  be  discussed. 

The  flexing  of  the  beam  adds  a  tremendous  complexity  to 
the  problem.  Referring  to  Fig  2,  it  can  be  seen  that  the 
bending  of  the  beam  will  change  the  overall  shape  of  the 
sh u t t 1 e -b e a m- a n t e nna  system  very  slightly.  An  assumption 
made  here  is  that  the  deformations  of  the  beam  is  small  along 
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Shat  tie 

Ij  (slug-ft2)  = 


mass  =  =  6366  .46 

905.443  0 

0  6,789.100 

-145.393  0 


slags 


-145,393 

0 

7,086.601 


Aatenni 

I4  (slug-ft2)  * 

aass  =  b^  ■  12.42  slags 

4.969  0  0 

0  4,969  0 

0  0  9,938 

Be  sb 

mass  *  bq  “  12.42  slags 

Roll  Beading: 

pA  “  0.09556  slags/ft 

El  -  4.0  x  107  lb-ft2 

£  =  0.003 

Pitch  Bending: 

pA  -  0.09556  slags/ft 

El  -  4.0  x  107  lb-ft2 

£  -  0.003 

Yaw  Torsion 

pj  -  0.9089  slug-ft 

GJ  -  4.0  x  107  lb-ft2 

£  -  0.003 

Proof-Mass 

Actaator  s : 

aass  *  *2  “  *3  “  0.3108  slags 
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its  entire  length,  and  thus  also  at  the  beaa-antenna  at¬ 


tachment  point. 

The  approach  to  setting  np  the  mathematical  model  in  this 
paper  will  be  to  assume  the  beam  and  shuttle  form  a  dynamical 
system  and  that  the  antenna  effect  is  to  subject  the  beam's 
free  end  to  forces  and  moments.  Since  the  antenna  is  assumed 
to  be  rigidly  attached  to  the  beam,  the  motion  of  the  beam  is 
totally  defined  by  the  rigid-body  motion  of  the  shuttle-beam 
system  plus  the  elastic  motion  of  the  beam  tip.  The  center 
of  mass  of  the  shuttle-beam-antenna  system  is  assumed  to  be 
unaffected  by  the  elastic  motion.  This  is  a  good  assumption, 
considering  that  the  mass  of  the  shuttle  is  over  500  times 
that  of  the  beam-antenna  combination  (see  Table  I). 

The  ot er al 1  mo t ion  of  the  beam-antenna  with  respect  to 
the  shuttle  is  to  be  controlled  by  force  and  moment  actuators 
located  on  the  antenna,  along  with  two  proof-mass  actuators 
(Fig  1),  located  at  positions  sn2  and  sn3  along  the  beam  to 
be  chosen  by  the  analyst.  The  optimum  location  for  these 
actuators  could  be  the  subject  for  an  entire  study.  For  this 
paper,  however,  the  actuators  will  be  assumed  to  be  located 
at  the  40-  and  80-foot  positions  along  the  beam  (see  Appendix 
A).  The  actuators  operate  by  moving  a  mass,  which  causes 
forces  in  the  x  and  y  directions  only.  There  is  no  torsional 
input  from  these  actuators.  There  is,  however,  a  moment 
created  about  the  system's  mass  center  which  will  tend  to 
rotate  the  shu 1 1 1  e -b e am- a n t e nna  system  out  of  its  desired 


attitude.  These  nonents,  although  small,  will  be  included  in 
this  study. 

Other  assumptions  made  are: 

1)  The  beam  does  not  appreciably  stretch,  meaning  that 
it  does  not  deform  in  the  z  direction. 

2)  Any  forces  in  the  z  direction  from  the  motion  of  the 
antenna  are  insignificant. 

3)  There  are  no  specific  forces  or  moments,  such  is 
those  due  to  meteor  collisions,  solar  rsdiation  pressure, 
gravity  torques,  or  magnetic  or  atmospheric  effects  modeled 
in  this  study.  These  forces,  if  they  exist,  will  be  small  in 
comparison  to  the  control  torques  available. 

The  next  chapter  will  develop  the  mathematical  model  of 
the  system,  which  in  its  general  form  would  be  nonlinear  and 
contain  both  partial  and  ordinary  differentials.  For  the 
study  at  hand  a  linear  discrete  model  was  produced  using  an 
assumed  modes  approach  and  by  considering  only  small  motions 
from  an  undeformed  equilibrium  position. 


This  section  will  sttack  the  mathematical  development  of 
the  system  model  in  fonr  parts.  The  first  vill  be  the 
choosing  of  functions  to  represent  the  flexing  of  the  beam. 
The  second  will  be  the  derivation  of  the  equations  of  motion 
of  the  entire  shut tl e-be am-antenna  system,  incorporating  the 
flexing  of  the  beam,  and  accounting  for  the  rigid  antenna 
through  forces  and  moments  acting  on  the  antenna  end  of  the 
beam.  The  third  will  be  the  development  of  the  equations  of 
motion  of  the  antenna  in  terms  of  the  displacements  and 
rotations  at  the  end  of  the  beam,  along  with  the  development 
of  the  forces  due  to  the  proof-mass  actuators.  The  fourth 
section  will  put  the  equations  of  motion  together  with  the 
generalized  forces  and  express  the  system  in  matrix  equation 
form  for  the  purposes  of  applying  a  control  law. 

Choosing  Proper  Functions 

A  discretization  approach  will  be  employed  in  this  paper 
in  choosing  functions  to  represent  the  motion  of  the  vib¬ 
rating  beam  with  respect  to  the  shuttle.  This  is  known  as 


an  assumed  modes  approximation.  The  more  modes  modeled,  the 


vibration  is  given  by  Meiroviteh  [1:208-209]  ss 

-32  EI(x)  32y(x,t)  +  f(x,t)  -  M(x)  32y(x,y)  (1) 

3  x2  3  x2  3t2 

If  the  beta  is  unifora,  this  is  reduced  to 

M  d2y(x.t)  +  El  d4y(x.t)  -  f(x.t)  (2) 

3t2  3x4 

where  y(x,t)  is  the  transverse  displacement,  x  is  a  coord¬ 
inate  along  the  beaa  length,  N  is  the  mass  per  unit  length, 
and  El  is  the  bending  stiffness.  For  free  vibration  the 
distributed  force  f(x,t)  *>  0.  Fox  the  case  of  a  cantilever 
beaa  the  associated  boundary  conditions  are 

y(0,t)  -  0 

3  y  (  x,  t )  I  “  0 
3x  lx«0 

EI32y( x, t )  -  0  (3) 

3x2  x=L 

EI3  3  y ( x , t )  -  0 

3x3  x  »L 

Denoting  the  running  length  of  the  bean  by  s,  the  aass  per 
unit  length  by  pA,  and  defining  the  two  orthogonal  coaponents 
of  bending  by  ur  and  up  (see  Fig.  3)  and  assuaing  the  beaa 


equations  of  motion  can  be  written: 


<0 


pA  d2ur(s,t)  +  El  34nr(s.t)  =  0 
at2  3s4 

pA  3 2Up ( s , t )  +  El  a4np(a,t)  =  0  (4) 

at2  as4 

subject  to  the  conditions 


ur(0, t)  *  up ( 0 , t )  =  0 
3ur<0,t)  *  dUp(O.t)  =  0 

a  s  a  s 


ei  a2ur(«, t) 

a.2 


-  ei  azuf>(s.  t) 


s*L 


a  s' 


s-L 


EI  d3ur(s.t) 


a  s 4 


s*L 


-  EI  d3up(s, t) 
ds3 


s=L 


(5) 


Since  the  two  deformations  ur  and  up  satisfy  identical 
relationships,  the  development  will  be  restricted  to  ur  and 
the  results  applied  to  both  ur  and  up.  Assume  that  ur  can  be 
written  as 


ur(s. t)  «  R(s)Dr(t)  (6) 

Substituting  eq  (6)  into  (4a)  and  dividing  by  the  product 
E( s ) Ur ( t )  yields 
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Ur  +  El  Riv  =  0 

Or  pA  R  (7) 

which  implies  that 

Riv  -  or4R  =0  (8 

where 

ar4  =  -pA  Dr  =  pA  ur2  (9 

El  Ur  El 

The  above  conditions  have  the  general  solntion 

R  =  Arsinars  +  Brcosars  +  Crsinhars  +  Drcoshars  (10 

Substituting  eq  (6)  into  the  boundary  conditions  yields 

R( 0)  *  R' (0)  =0  (11 

R'  '(L)  =  R'  "(L)  =  0 

Solving  these  four  equations  simultaneously,  the  results  are 
three  important  constraints  which  drive  the  rest  of  the  dev- 
e 1 opment : 
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Ar  =  -Cr 

-Dr  =  -Ar(sinurL  +  sinhazL) 


(coso.L  +  cosha.L) 


( 12 


The  lest  of  eqs  (12)  has  an  infinite  number  of  solutions. 
These  solutions  will  be  denoted  by  ari  an<i  the  associated 
amplitudes  by  Ari,  Bri.  Cri.  and  Dri-  The  mode  shapes  are 
then  given  by 

Rri  =  Ari  [<sinori8  ~  si«*“ris) 

+  (sinor£L  +  s  i  nhar  jL)  (  co  shor  £  s  -  cosar£s)J  (13) 

(cosariL  +  coshar^L) 


The  amplitude  Ari  is  arbitrary.  For  convenience,  the 
magnitude  Ar^  will  be  chosen  such  that 

f  pA  R^R ^ds  =  1  (14) 

J  0 

Expanding  (14)  and  integrating  each  individual  term  dir¬ 
ectly  takes  a  great  deal  of  bookkeeping  and  substitution  (see 
Appendix  B),  but  the  solution  reduces  to: 


kri  =  /  1  N1^2  <c°*ariL  +  co 

\(pAL)/  (sinar£L  +  si 


co  shor  ^L) 


sinhdj.  ^L) 


B 


r  l 


f— ) 

VpAL)/ 


1/2 


(15) 


Since  the  mode  shapes  have  been  normalized,  they  possess  the 
convenient  property 


f 


pA  R  jRj  ds  =  6  £  j 


(16) 


where  6  ,  { 


is  the  Kronecker  delta. 


The  proceeding  derivation  applies  to  a  cantilever  beam 
in  bending.  The  beam  element  in  the  model  of  Fig  (3)  is 
assumed  to  behave  as  a  cantilever  in  bending  in  each  of  two 
orthogonal  directions  (roll  and  pitch),  and  in  addition  to 
undergo  torsional  deformation.  The  normal  modes  given  by  eq 
(13)  may  also  be  used  for  pitch  bending  if  the  subscripts  r 
are  replaced  with  p  throughout.  The  roll  motion  is  therefore 
given  by: 


X 

‘r  *  *i< 


s)Uri(t) 


and  the  pitch  motion  as 


up  =  S  Pi<* 


)  Up  i  ( t ) 


where  nr  is  the  number  of  roll  modes  and  np  is  the  number  of 
pitch  modes. 

The  torsional  motion  about  the  z  (yaw)  axis,  however, 
must  be  treated  separately.  The  differential  equation  for 
the  torsional  motion  on  a  beam  is  given  by 


pJ32uy  -  GJd2ny  =  0 


where  J  is  the  polar  area  moment  of  inertia.  The  associated 
boundary  conditions  are  given  by 


uy ( 0, t )  =  0 

u'y(L. t)  =  0 
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Proceeding  as  in  the  case  of  bending  vibration,  the  torsional 
mode  shapes  are  given  by 


where 


Yi  =  CyisinPyis 


0yiL  =  (2n-l)n/2 


n  —  1,2,3,. 


As  before,  it  is  convenient  that  the  mode  shapes  be  ortho¬ 
gonal,  so  the  condition  that  mast  be  satisfied  is 

e  L 


(  PJ  TjY, 
j  0 


ds  -  1 


Substituting  (21)  into  (23)  results  in 
r  L 

l  pj  Cyi2sin2py.s  ds  =  1 
*  0 

which,  after  integrating  and  applying  the  limits,  gives 


\pJ(20y.L  -  sin2pyiL)/  \ 

pJL  / 

The 

sine  term  goes  to  zero  because  2*beta*L 

will 

result 

in  an 

integer  multiple  of  pi  for  all  betas. 

As 

before,  these  functions  demonstrate  the 

property 

rL 

J0  pJ  YiYi  ds  -  6ij 

(26) 

The 

torsional  motion  may  now  be  described  by  the 

rel ationship 

y ' 

®y(S,  t)  =  2~!  Yi<«)U_i(t) 

(27) 

aiL 

(radians) 

Ai 

Bi 

n2 

2 

( r ad / sec) 

1.875102 

.2082777241 

-  .2837201974 

18.11807 

4 .694091 

.2889597487 

-.2837201974 

711.56816 

7 .854757 

.2835001714 

-  .2837201974 

5578.80978 

10.995541 

.2837297171 

-  .2837201974 

21422 .82836 

14.137168 

.2837197  860 

-.2837201974 

58541  .00381 

17 .278759 

.2837202151 

-  .2837201974 

130635.34107 

20.420352 

.2837201966 

-.2837201974 

254837 .51330 

23 .561945 

.2837201974 

- .2837201974 

451705 .09099 

26 .703538 

.2837201974 

-.2837201974 

745221 .94382 

29.845130 

.2837201974 

- .2837201974 

1162798 .20573 

32.986723 

.2837201974 

-  .2837201974 

1735270 .27761 

36 .128315 

.2837201974 

-  .2837201974 

2496900 .82710 

39 .269908 

.2837201974 

-.2  837201974 

3485378.78863 

42.411501 

.2837201974 

-.2837201974 

4741819 .36334 

where  ny  is  the  number  of  yaw  modes,  Uyi(t)  are  time  depen¬ 
dent  modal  amplitudes  and  the  Y^(s)  are  given  by  eq  (21). 

A  computer  program  was  written  to  calculate  the  roots  of 
eq  (12c).  Using  the  first  fourteen  roots,  the  first  fourteen 
roll  and  pitch  functions  were  evaluated,  solving  for  each 
and  the  (constant)  B^s.  The  u^'s  were  also  computed,  and 
the  data  is  tabulated  in  Table  II.  The  betas,  being  multi- 
pies  of  pi,  were  found  directly  and  the  Cyj's  and  wy  ^  '  s  were 
calculated  for  the  yaw  torsion  functions.  This  data  is 
presented  in  Table  III. 

In  the  next  section,  certain  integral  relationships  in- 


Table  III 

Angles.  Coefficients,  end  Squared  Frequencies  for 
Taw  Torsion  Equation 


PiL 

(radians ) 

Ci 

-i2 

(rad/sec)^ 

1 .570796327 

.1301023886 

6425.3522 

4 .712388980 

.1301023886 

57828 .1697 

7 .853981634 

.1301023  886 

160633 .8047 

10.9955742 9 

.1301023886 

314842.2572 

14.13716694 

.1301023866 

520453 .5273 

17 .27875959 

.1301023  866 

777467.6148 

20  .42035225 

.1301023866 

1085884 .520 

23 .56194490 

.1301023  866 

1445704.242 

26  .70353756 

.1301023  866 

1856926  .782 

29.84513021 

.1301023  866 

2319552 .140 

32.98672286 

.1301023866 

283  3  580  .315 

36 .12831552 

.1301023866 

3399011  .308 

39.26990817 

.130102- 866 

4015845.118 

42 .41150082 

.1301023866 

4684081 .745 

volving  the  continuous  coordinates  ur,  np,  and  uy  will  be 
needed.  Specifically,  the  kinetic  and  potential  energies 
associated  with  these  variables  are  of  interest.  These  ener¬ 
gies  involve  the  following  six  integrals: 


functions  for  s  fixed-free  uniform  rod  in  both  bending  end 
torsion.  These  functions  will  be  nsed  es  assumed  modes  for 
the  beam  part  of  the  shut 1 1  e-beam- ant e nna  system  in  the  next 
section.  The  properties  of  these  functions  given  in  eqs  (29) 
and  (30)  will  prove  useful  in  this  development. 

System  EqyiUpaj 

The  system's  equations  of  motion  will  be  developed  using 
a  combination  of  Lagrange's  equations  for  the  elastic  motion 
of  the  beam  and  Euler's  moment  equations  for  the  overall 
motion  of  the  shuttle-beam  combination.  To  that  end  the 
kinetic  and  potential  energies  for  the  shuttle-beam  combina¬ 
tion  will  be  developed.  The  effect  of  the  antenna  will  be 


taken  into  account  through  the  generalized  forces  due  to  the 


forces  end  noments  present  at  the  bean-antenna  attachment 


point.  The  b  i 


>d  to  be  rigidly  attached  to  the 


shuttle  at  the  center  of  si  a  s  s  of  the  system.  It  is  further 


long  axis  and  to  undergo  torsional  notion  around  its  long 


axis. 


The  location  of  a  general  point  in  the  shuttle-beam 


system  relative  to  the  system's  mass  center  is  given  by 


R  =  r  +  n 


where  r  =  xx  +  yy  +  zz  locates  a  generic  point  and  u  denotes 
the  elastic  displacement  of  the  mass  particle  at  r.  Defining 
R,  and  Rb  as  position  vectors  of  points  in  the  shuttle  and 
beam  respectively,  the  equation  becomes 


Rj  *  xx  +  yy  +  zz 

Rb  *  (  u  j.  (  z.  t)  +  x )  x  +  ( Up  (  z .  t )  +  y)  y  +  zz 


The  system  kinetic  energy  is  thus  given  by 


T  =  1/2  «tVc2  +  1/2  j 


.  1  .  1 

R,  Rs  dm. 


f  .  I  .  I 

+  1/2  /  R„  -Rb  dmb 


where  m^  is  the  total  mass,  Vc  is  the  velocity  of  the  mass 


tives.  Denoting  the  angular  velocity  of  the  body-fixed  shut¬ 
tle  axes  by  u  and  that  of  an  eleaent  of  the  beam  by  ,  the 
kinetic  energy  can  be  rewritten  as 


1/2  “tvc2 


+  1/2  J  ( u  x  Rs)*(w  x 


Rs)d-s 


(34) 


f  ,  B 

+  1/2  /  (Rb  + 


B 


“b  x  Rb^'^®b  +  "b  x  *b^d"b 


where  the  superscript  B  denotes  derivatives  seen  by  an 
observer  fixed  in  the  shuttle  body  axis,  so  R  8  ®  -  0  since 
the  shuttle  is  rigid.  Furthermore  it  should  be  noted  that 

.  B 

Rb  =ur(x,t)x  +  ip(z,t)y  (35) 


and  that  to  and  wb  are  related  by 

“j,  =  u  +  «^/s  *  •  +  iy(**t)x 


(36) 


where  uy(z,t)  is  the  torsional  angular  displaceaent  of  the 
beam  eleaent  at  position  z.  The  kinetic  energy  is  expanded 
and  expressed  in  matrix  notation  as 


) 


+ 


+ 


where  Is  is  the  sonent  of  inertia  aatrix  for  the  shuttle. 
The  lattices  u  and  /  s  *n  e<*  (37)  are  foried  froa  the 

coaponents  of  the  u  and  ^  s  vectors  in  the  form 


[s] 


(38) 


The  kinetic  energy  in  eq  (37)  contains  teras  through  order 
fonr  in  the  variables  ur,  up,  «i>y,  and  uy  and  their  rates. 
Since  it  is  desirable  to  provide  a  linear  set  of  equations 
for  the  aotion  of  the  system  in  the  neighborhood  of  an  equi¬ 
librium  where  all  of  these  variables  and  their  rates  are 
small,  T  will  be  reduced  to  include  only  those  terms  of  order 
two  or  less  in  these  variables.  This  yields  the  kinetic 
energy  as 

T  ■  Ill  |»|T  [i,  +  b]  {.} 
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39 


1*M 

tB 

( 


Expansion  of  eq  (39)  yields  numerous  terms  including  integ¬ 
rals  involving  the  variables  x  and  y  linearly  and  the  combi¬ 
nations  xy,  xz,  and  yz.  Since  the  elastic  deformations  uc. 
Up,  and  Uy  spatially  depend  only  on  z,  these  integrals  are 
all  zero  from  symmetry  considerations.  Setting  these  terms 
equal  tozeroyields  the  kinetic  energy  in  the  form 


T 


1/2 |«}T 


jui  +  1/2 


A  0  0 
0  A  0 
0  0  J 


ur  \ 


z 

0 


0 

0 


0  x  2  +  y  2  J 


j  u»  |  dz 


dz 


(40) 


Substituting  the  assumed  modes  from  the  previous  section 
results  in 


Three  blocks  of  terms  in  the  Sz  matrix  warrant  closer 


attention. 


The  lover  right  corner  block  is  x2  +  y2  and. 


after  aatrix  ■ nl t ipl i ca t ion  and  integration,  will  result  in 


Sz(3,k)  =  ^  2W3J  Uy £ dz 


(43) 


where 


i  *  1,  2  ,  3,  ..., 


k  *  nr  +  np  +  1 


These  terns  are  considered  insignificant  for  this  analysis, 
since  a  typical  value  of  J  will  nake  the  terns  snail  conpared 
to  the  others  in  the  natriz.  Therefore  these  terns  will  be 
assnned  to  be  zero. 

The  other  non-zero  terns  in  the  natrix  are 


Sz  ( 1 ,  j  )  =  pAzPj(z)dz 


(44) 


and 


Sz  (  2  ,  i )  =  -  Sz  ( 1 ,  j  ) 


(45) 


where 

i  —  1,  2,  3,  ...,  Hji  j  =  n£  +  1,  nr  +  2,  ...,  n r  +  Up 

Substituting  the  eigenfunction  expressions  into  (44)  yields 

fL 

Sz(l.j)  =  -pA  j  [zApj s inapj z dz  +  zBp j co sap j z dz  (46) 

-  zAp j s i nhap j z dz  -  zBp j co shap j z] dz 
Integrating  each  of  the  four  terns  separately,  and  taking 


advantage  of  the  fact  that 


Bpj(cosopjL  +  co*hapjL)  *  -Apj(sinapjL  +  sinhapjL)  (47) 

and 


Bpj  =  -1/  (pAL)1^ 


(48) 


eq  ( 46 )  be  cone  a 


Sx(l.j)  =  2 ( p AL^ ) !/ ^ 
(°PjL)2 


(49) 


The  Sz(2,i)  ter*  will  differ  by  an  algebriac  sign,  and  the 
snbscripts  will  becoie  ri  instead  of  pj .  The  potential 
energy  can  be  expressed  as 


(  50) 


Using  the  assaned  node  expressions,  this  bee ones 


V 


1/  2 


So,  forming  the  Lagrangian  and  taking  derivatives. 


(51) 


L  =  T  -  V 


3  L 

a|u} 


(  52) 
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Forming  Lagrange  $  equations : 


d 

dt 


d  L 


aDi 


(53) 


where  Q  ^  are  the  generalized  forces  to  he  developed  in  the 
next  section.  Substituting  in  the  terms  just  developed.  La¬ 
grange's  equations  for  this  problem  become 


(54) 


This  equation  assumes  no  damping  in  the  beam.  The  damping 
will  be  figured  in  as  one  of  the  last  steps  before  applying 
the  control  law. 

Eq  (54)  is  but  one  of  two  equations  which  can  be  obtained 
from  the  Lagrangian  method.  If  derivatives  are  taken  with 
respect  to  <■>,  the  result  is 


3T 


d  <t> 


Applying  Lagrange's  method,  the  result  is 


(55) 


(56) 


where  N  is  the  sum  of  all  moments  on  the  system.  The  second 
term  on  the  left  side  of  the  equation  can  be  ignored  since  it 


is  of  higher  order.  The  remaining  terms  give  the  equation 

d  J  9T  =  |  M  j>  (: 

dt  3<i) 


which  are  generally  referred  to  as  Lagrange's  equations  in 
q ua s i- c o or d i na t e s.  Again,  substituting  in  the  previously 
derived  terms,  the  equation  becomes 


]  W  +  [s*]  M  ■  W 


This  equation  along  with  eq  (54)  will  serve  as  the  two  equa¬ 
tions  of  motion  for  the  entire  shut tl e-beam- antenna  system. 
Eq  (54)  is  the  rigid  body  equation  of  motion  with  a  modifica¬ 
tion  taking  into  account  the  flexing  of  the  beam.  Eq  (58)  is 
Euler's  moment  equation  with  a  coupling  term,  again  to  ac¬ 
count  for  the  beam's  flexing. 

Now  that  the  functions  have  been  chosen  and  the  two  equa¬ 
tions  of  motion  derived,  it  is  time  to  turn  to  the  develop¬ 
ment  of  the  forces  and  moments  on  the  cantilever  beam.  This 
next  section  will  focus  on  the  term  by  term  derivation  of 
these  forces  and  moments. 


Now  that  the  unforced,  undamped  equations  for  the  roll, 
pitch,  and  yaw  bending  have  been  obtained,  the  next  step  is 
to  develop  the  generalized  forces  and  the  moments  N  of 
eq  s  ( 5  4 )  and  (58). 

The  generalized  forces  are  determined  from  [4]  as 
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g 


a 


j  =  2  aui  k=2  aui 

_  • 
where  the  f,  are  the  applied  forces,  the  r,  is  the  velocity 
J  J 

of  the  point  of  application  of  the  force,  and  the  D  ^  are  the 
modal  velocities  from  eq  (42),  where  i  runs  from  1  to  nr  + 

+  ny.  The  g^  are  the  applied  torques  and  is  the  angular 
velocity  of  the  element  at  which  the  torque  is  applied.  For 
the  problem  at  hand  the  only  applied  forces  of  interest  are 
those  shown  in  Fig  4  due  to  the  actuators  and  the  forces  and 
moments  due  to  the  antenna.  From  Fig  4  it  can  be  seen  that 
the  forces  and  moments  are  applied  at  three  specific  loca¬ 
tions.  Specifically,  the  forces  are  applied  at  points  sn2 
(which  is  40),  sn3  (which  is  80),  and  at  the  beam-antenna 
attachment  point  sn4  (which  is  130).  Notice  that  the  index  j 
ranges  from  2  to  4.  This  is  because  all  forces  at  snl  (the 
shuttle-beam  attachment  point)  are  zero  due  to  the  cantilever 
model. 

Moments  are  applied  only  at  sn4.  As  before,  all  moments 
at  snl  are  zero,  so  the  index  k  ranges  from  2  to  4.  Each  of 
the  forces  and  moments  from  eq  (59)  must  now  be  individually 
identified. 

The  proof-mass  actuators  are  designed  to  apply  a  force  in 
the  X  and  Y  directions  only.  The  X-direction  forces  are: 


I 


fr,2  =  “2  d2nc  +  Fr, 2 

3t2  s= 40 

f r, 3  =  ®3  ^  2pr 

at2 

where  the  first  term  in  each  equation  is  the  mass  of  the 
actuator  multiplied  by  the  acceleration  of  the  be  am- actua tor 
attachment  point,  and  the  second  term  is  the  force  metered 
out  by  the  actuator  itself.  The  subscript  2  or  3  denotes  the 
actuator  at  position  sn2  or  position  sn3  respectively.  The 
T-direction  forces  have  an  identical  form: 

f p , 2  =  “2  a2°p 

at2 

fp,3  =  “3  *2up  +  Fp , 3  (61) 

at2  s=  80 

The  last  forces  and  moments  to  contend  with  are  the  forces 
and  moments  on  the  beam-antenna  attachment  point.  The  NASA 
paper  [2]  gives  expressions  for  these  forces  which  appear  to 
be  in  error.  Therefore  the  expressions  for  the  forces  and 
moments  are  derived  as  follows. 

Looking  at  the  free-body  diagram  of  the  antenna  (Fig  4). 
the  forces  will  be  developed  from  the  force  equation.  The 
position,  velocity,  and  acceleration  of  the  center  of  mass  of 
the  antenna  with  respect  to  the  attachment  point  p  are: 


+  Pn  , 

p  ,  2 

s=40 
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r c  =  18 .75*4  -  32 .5y4 


vc  =  drc  +  “4  x  rc  =  “4  x 


*  c  =  d(u>4  i  r  c )  +  u4  x  (w4  x  rc)  =  £>4  x 


These  equations  make  use  of  the  fact  that  the  antenna  is 
modeled  as  a  rigid  body,  so  the  time  rate  of  change  of  the 
position  vector  in  this  frame  is  zero.  All  higher  order 
terms  have  been  neglected.  Now,  using  F=MA,  and  denoting 
the  acceleration  of  the  center  of  mass  of  the  antenna  in  this 
frame  as  Ic  : 


f 4  =  ®4ac  “  m4(ip  +  ac/p>  (63) 

where  ip  is  the  acceleration  of  the  attachment  point  and  ac/p 
is  the  acceleration  of  the  mass  center  with  respect  to  point 
p.  The  acceleration  of  point  p  is,  in  vector  form: 


ap  ■  a2*rx  +  a2upy 


The  acceleration  of  the  mass  center  wrt  point  p  takes  a 
little  more  development.  Since  the  antenna  is  rigidly 
attached  to  the  beam,  <i>,  xs  given  by 


d  sdt  s=13 


“4  * 


•s 

+  (0  = 

“4b 

9  s  9 1 

s=  1 3  0 

8ay 

“4c 

9 1 

s=  1 3  0 

+  a) 


(65) 


where  <o  is  the  rotation  of  the  entire  sh  n  1 1 1  e- b  e  a  m- an  t  e  nn  a 
system.  This  effect  was  addressed  in  the  rigid  body  equation 
derivation  and  will  not  appear  in  the  ant e nna- f ixe d  reference 
frame.  Therefore,  the  time-derivative  of  is  found  to  be: 


u 

0 

CO 

1 _ 

9  sdt2 

s=  1 3  0 

»S 

3s9t2 

s=130 

s2®y 

at2 

s=  1 3  0 

(66) 


which  can  be  written  in  vector  form  as: 


« . 


=  a 


£  +  9  ^  n 


A  .  ,2 

y  +  a^u 


9  s  9 1  • 


3  s  9  t2 


9  t  ‘ 


i=130 


(67) 


Before  taking  the  cross  product  of  the  and  the  posi¬ 
tion  vector,  the  vectors  must  be  in  compatable  frames.  It 
can  be  shown  that  for  very  small  angles,  the  frame  and 

the  xyz  frame  are  essentially  equal.  This  derivation  will 


A 


A 


1 


< 


« 


A 


obtain  the  same  resalt,  bat  will  go  through  the  cros s-prodact 
step  by  step,  showing  what  assumptions  are  necessary. 

Referring  to  Fig.  (5),  three  possible  rotations  oat  of  the 
x4Y4z4  *re  shown,  along  with  their  respective  rotation  mat¬ 
rices.  Using  the  small  angle  assumption  that  the  cosine  of 
an  angle  is  approximately  equal  to  1  and  the  sine  is  approxi¬ 
mately  equal  to  the  angle  itself,  these  three  matrices  become 


1  -e. 


R1  -  ®4  z 


R->  = 


o  e4y 


r3  = 


Multiplying  these  together  and  ignoring  any  non-linear  terms, 
the  rotation  matrix  from  the  x4y4z4  frame  into  the  xyz  frame 


1  "e4z  ®4y 


-04y  e4x 
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where  the  Rx/4  denotes  rotation  of  the  147424  frame  with 
respect  to  the  xyz  frame.  Therefore,  the  position  vector 
becomes,  in  colnmn  vector  form 


18.7  5  +  32.5  0, 


ic  =  18.7504z  -  32.5 


-18.7504y  -  3  2 . 5  O4  x 


Taking  the  cross  product  using  matrices,  the  a>4  vector  is 
written  in  'tilda'  form  and  the  cross  product  becomes 


0*  • 

A  „  0). 


*4  c  “4b 


“4c  0  "“4a 


18.75  +  3  2 . 5  O4  z 


1 8 . 7  50x  _  -  32.5  (70) 


~“4b  “4a  ®  \-l8.75©4y  -  32.5©4X 


Therefore,  the  acceleration  of  the  mass  center  wrt  point  p 
is  given  by  (ignoring  non-linear  terms): 


ac/p  =  3  2 . 5  «i>4  cx  +  18.75«4cy  +  (-18.75w4b  -  32.5<L4a)z 
=  3  2.53  2Uy.x  +  I8.7532uyy  -  /32.533ur  +  18.7533u\£  (71) 


3  s3t2 


3  s  3 14 


The  z  component  will  be  ignored,  since  the  assumption  is 
that  the  beam  is  not  stretched  in  the  z-direction.  So, 
putting  all  of  the  terms  together  in  the  force  equation, 

f,  -  m^  f/  3  2u.  +  3  2.53  2u„\  x  +/d2u„  +  I8.7532u„\y1  (72) 


Breaking  these  forces  into  their  respective  components 
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r4 
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p4 
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s=  1 3  0 
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s=130 


32.5m4  a2uy 
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These  are  all  of  the  forces  to  be  modeled  in  this  stndy.  The 
expressions  for  the  moments  will  now  be  developed. 

The  only  moments  to  derive  are  those  about  the  beam- 
antenna  attachment  point  (see  Fig.  4).  These  are  obtained 
from  Euler's  moment  equations  and  are 


gr ,  4 
gP.  4 
gy .  4 


(73) 


The  middle  term  will  be  neglected  because  it  is  of  second 
order.  The  u>4  was  given  by  eq  (6  5)  as 


“4 


a2u. 


a  sat 


a2*. 
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s= 1 3  o 
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s  =  1 3  0 


(65) 


Taking  the  time  derivative  of  eq  (65)  and  substituting  in  the 


mode  approximations,  the  result  is 


UrR' ( 13  0  V 
“4  =  {  DpP* (130) 


0y Y (130) 


From  Table  I,  1-4  is  given  as 


(74) 


4969  0  0 

0  4969  0 

0  0  993  8 


So  the  moments  become 


(  gr ,  4 

I  ®p,  4 

'  8y ,  4 

4969  0  0 

0  4969  0 
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UrR' (130) 

UpP '(130))  + 
\  DyY(l30) 
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(75) 


which  reduces  to 


4  96  9UrR* (130)  +  M4i j 

4969UpP' (130)  +  K4yl  (76) 

99380yY(130)  +  M4  z  ] 

Now  that  each  force  and  moment  has  been  identified,  each 
Qj  must  be  developed.  A  sample  derivation  of  follows. 

Every  other  Qj  ranging  from  Q2  through  Q(  nr  +  ^  +  ny )  wil1  have 
a  similar  development. 

From  eq  (59),  the  will  be  given  by 


The  r2  term  is  derived  from 


t2  =  (x  +  ur)x  +  (y  +  u  )y  ~  sn2z 


r  P 

r2  =  5  x  r2  +  ^DiR^dOJx  +  Pj  (40)y 


So  the  partial  derivative  with  respect  to  P,  is 


dr2  =  1^40)1 


S imil arly. 


di3  =  R1(80)x 


3r4  =  Rj(130)x 


The  angular  velocity  u>4  was  given  by  eq  (65).  Substituting 
the  assumed  modes  yields 


_ *  y 

5a  =  X]p4R4  '  (130)x  +  TVp4.<130» 


tac“  01  tlie  or  +  np+ny  generalized  forces  can  be  derived  into 
similar  expressions.  Since  for  this  model  14  modes  for  each 
of  the  roll,  pitch,  and  yaw  motions  are  assumed,  there  are  42 
Q  £ 1  a  in  all.  The  first  14  will  have  terms  associated  with 
the  second  derivative  of  the  roll  amplitudes,  terms  asso¬ 
ciated  with  the  second  derivative  of  the  yaw  amplitudes,  and 
terms  associated  with  the  proof-mass  actuator  forces  in  the 
x-direction,  as  shown  by  eq  (85).  Expanding  these  terms  in 
matrix  form  results  in  three  matrices.  For  example,  using 
just  the  14  roll  terms,  the  matrix  associated  with  the  second 
derivative  of  the  roll  amplitudes  is 


trrl.l 
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where 


trri>j  =  -  m2 R i ( 40 ) R j ( 40 )  -  m3  R£  (  80  )  Rj  (  80  ) 

-  m4Ri(130)Rj (130)  -  4  969 R  £  '  <  13  0  )  R j  ( 1 3  0 ) 

The  matrix  associated  with  the  second  time  derivative  of  the 
yaw  amplitudes  looks  like 
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where 

tryi#  j  =  -32.5B4Ri(  130)  Yj  (  130) 
and  the  matrix  of  forces  and  moments  is 

-Fr,2Rl(40)  -  Fr>3R1(80)  -  M4iR1'(130) 

~Fr,2R2(40)  “  Fr,3R2(80)  "  M4xR2'{130) 

Trc  =  -Fr>2R3(40)  -  Fr>3R3(80)  -  M4xR3'(130)  (88) 

•  •  • 

«  •  • 

•  •  • 

-Fr  2R14(80)  -  Fr3R14(80)  -  M4xR14'(130) 

The  pitch  equation  will  have  a  right  side  in  an  identical 
form  with  the  only  change  being  that  any  roll  functions  are 
replaced  by  pitch  functions.  The  yaw  equation  simplifies 


ciated  with  the  second  time  derivative  of  the  yaw  amplitudes 
and  a  matrix  of  force  and  moment  terms.  There  is  no  coupling 


in  the  yaw  equation  with  either  pitch  or  roll,  which  is 
expected  for  a  fixed-free  model.  The  yaw  matrices  are 
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where 


tyyi.j  =  -7.086,601Yi(0)Yj  (  0)  -  993  8Y£  ( 13  0  )  Yj  (  13  0  ) 

and 


-M4zYl( 130) 
-M4zY2(130) 

-M4zY3(130) 

-M4zYi4(13°) 


(  90) 


These  matrices  come  from  the  last  14  Qz's,  that  is,  Q2  9 
through  Q42  . 


Final  Model 

Before  tying  together  the  equations  of  motion  incorpor¬ 
ating  the  generalized  forces,  eq  (54)  will  be  slightly  modi¬ 
fied  to  represent  the  angular  rate  of  the  shuttle-beis- 
antenna  system  and  its  time  derivative  in  terms  of  the  corre¬ 


sponding  angular  displacements  and  their  time  derivatives. 
The  middle  term  on  the  left  side  of  eq  (54)  can  be  written  as 


where  ©j,  62*  and  63  are  the  angles  through  which  the  shuttle 
system  rotates.  Similarly,  the  first  term  on  the  left  side 
of  eq  (58)  can  be  rewritten  as 


The  left  side  of  eq  (54)  is  in  a  very  simple  form,  thanks 
to  the  judicious  choice  of  assumed  modes.  The  right  side, 
however,  contains  terms  in  the  generalized  forces  which  will 
combine  with  terms  on  the  left  side  and  complicate  them 
somewhat.  Specifically,  each  of  the  generalized  forces  con¬ 
tains  terms  which  will  combine  with  the  identity  matrix  (the 
mass  matrix)  on  the  left  side  of  the  equation.  If  eq  (54)  is 
rewritten  to  incorporate  the  generalized  forces,  it  can  be 
expressed  as 
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where  each  single  term  in  the  matrices  accociated  with  the 
second  derivative  of  the  amplitudes  represents  a  14  x  14 
block,  and  each  I  represents  a  14  x  14  identity  matrix.  The 


first  matrix  on  the  right  side  can  be  combined  with  the  first 
matrix  on  the  left  side.  The  resulting  equation  can  be 
written  as 


where 
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(94) 


A  new  vector  must  now  be  defined  to  augment  the  U  vector. 
It  is  formed  by  incorporating  the  O' s  and  can  be  written: 


X  =  |e  U,  Dp  Uy}T  (95) 

This  is  vector  of  45  elements  with  the  first  3  elements  being 
the  angles  through  which  the  shut t 1 e-b e am- a nt e nna  system 
rotates.  Eq s  (94)  and  (58)  can  now  be  added,  and  their  sum 
e  xpr e  s  se  d  as 
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remembering  that  Sz  is  not  a  square  matrix  (in  fact,  it  is  a 


*"  '  ■V - T — 3— 


v*  v»  v*  v»  •.  .  ■.  x  . 


3  by  42  for  a  14-mode  approximation).  The  new  mass  matrix  in 
this  equation  will  now  be  called  M  (no  subscript),  and  the 
force  and  moment  matrix  on  the  right  side  will  be  called  F. 
The  stiffness  matrix,  which  now  has  zeros  as  the  first  three 
diagonal  elements,  will  continue  to  be  called  co^. 

One  final  operation  will  be  performed  on  this  equation 
for  convenience.  It  is  desirable  for  the  mass  matrix  to  be 
the  identity  matrix  as  it  was  prior  to  the  incorporation  of 
the  generalized  forces  so  that  the  state-space  form  will  be 
easier  with  which  to  work.  It  is  possible  to  accomplish  this 
and  still  keep  the  stiffness  matrix  in  diagonal  form  by  using 
some  of  the  properties  of  matrix  manipulation.  If  an  eigen¬ 
value  analysis  is  done  on  the  unforced  (homogeneous)  form  of 
eq  (96),  the  resulting  eigenvectors  can  be  put,  column  by 
column,  into  a  transfer  matrix  {5:182-1861.  The  vector  x  can 
be  expressed  as 

i  =  0  ii  (97) 
so  eq  (96)  can  be  rewritten  as 


MMN  ♦  MIMM  -  H 


(98) 


Now,  pr e-mul t ip ly ing  eq  (98)  by 
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(  99) 


the  results  will  be 
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(100) 
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where 


K  =  M  -1  0  -1  N(u 2  0 


-i  U  -i 


( 101) 


The  mass  matrix  has  thus  been  reduced  to  the  identity 
matrix,  and  the  new  stiffness  matrix,  K,  keeps  its  diagonal 
form.  The  terms  on  the  diagonal  of  K  also  happen  to  be  the 
eigenvalues  of  the  unforced  system.  The  term  on  the  right 
side  of  the  equation  will  now  be  called  EF,  and  will  figure 
prominently  in  the  control  portion  in  the  next  chapter. 

Finally,  the  damping  of  the  beam  must  be  taken  into 
account.  The  damping  matrix  can  be  defined  as 


[D]  -  [2  5-] 


(102) 


where  the  £  was  defined  in  the  reference  [2]  as  .003  for 
roll,  pitch,  and  yaw  motion.  This  leaves  the  final  mathema¬ 
tical  representation  of  this  system  as 


'■]{*}♦  [e]M 


(103) 


Each  matrix  on  the  left  side  is  a  45  x  45  diagonal  matrix, 
and  the  matrix  on  the  right  side  is  a  45  x  1  matrix,  which 
will  be  modified  in  the  next  chapter. 

The  mathematical  model  of  the  shuttle-beam-antenna  system 
has  now  been  derived,  resulting  in  a  single  matrix  equation. 
The  next  chapter  will  apply  apply  linear  control  theory  to  it 
to  investigate  how  the  system  might  be  controlled. 
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The  previous  chapters  have  taken  a  complicated  physical 
system  and  reduced  it  into  a  mathematical  model.  This  model 
has  simplified  the  system  somewhat  by  ignoring  numerous  non¬ 
linear  terms.  The  final  mathematical  expression  consists  of 
a  fairly  large  matrix  equation  which  will  now  be  used  to 
develop  a  control  law. 

The  bulk  of  the  control  work  done  in  this  paper  is  a 
direct  result  of  in-depth  analyses  done  by  Janiszewski  [6] 
and  Aldridge  [7].  The  computer  programs  generated  by  this 
work  were  modified  for  use  on  this  model,  and  other  than 
observing  results,  no  attempt  was  made  to  further  the  study 
of  the  control  techniques  employed.  A  brief  outline  of  the 
theory  will  be  presented  in  this  chapter.  However,  anyone 
wishing  detailed  study  of  this  control  method  should  refer  to 
the  works  referenced  throughout  the  chapter. 

Since  it  would  be  impossible  to  control  all  of  the  modes 
of  a  vibrating  structure,  a  method  must  be  found  to  control 
the  most  excitable  modes  while  being  careful  not  to  drive  any 
uncontrolled  modes  unstable.  Janiszewski  [6]  shows  how  this 
is  done  by  dividing  the  modes  into  three  categories:  control¬ 
led,  suppressed,  and  residual.  The  controlled  modes  will  be 
actively  controlled.  The  suppressed  modes  will  not  be  ac¬ 
tively  controlled,  but  care  will  be  taken  to  avoid  exciting 


them  while  working  on  the  controlled  modes.  The  residual 
modes  will  not  be  controlled  or  suppressed,  with  the  as¬ 
sumption  that  their  frequencies  are  either  too  high  to  excite 
significantly  or,  if  excited,  will  dampen  out  through  the 
natural  damping  of  the  beam. 

The  first  step  towards  controlling  the  system  is  to 
truncate  the  mathematical  model  to  a  reasonable  number  of 
modes  with  which  to  work.  The  ACOSS  program  developed  and 
modified  by  Janiszewski  and  Aldridge  could  be  easily  modified 
for  this  model  and  handled  twelve  modes.  Since  the  modes 
with  the  lowest  frequencies  tend  to  be  the  most  excitable, 
the  model's  twelve  lowest  modes  were  used.  To  get  just  those 
modes,  eq  (103)  was  put  in  the  proper  form  by  judicious  use 
of  the  0  matrix.  The  eigenvectors  of  the  eigenvalue  problem 
solution,  which  are  the  columns  of  ©  ,  were  ordered  so  that 
the  eigenvalues  which  appear  on  the  diagonal  of  the  K  matrix 
were  in  ascending  order.  The  three  lowest  frequencies  (the 
square  root  of  the  eigenvalues),  which  are  zero,  correspond 
to  the  rigid  body  modes  and  must  necessarily  be  controlled  to 
keep  the  shuttle  in  the  proper  attitude.  The  next  nine  modes 
are  the  lowest  of  the  roll,  pitch,  and  yaw  beam  modes,  or¬ 
dered  from  lowest  to  high  est  without  regard  to  which  axis 
they  correspond.  It  is  these  modes  which  are  the  most  ex¬ 
citable  and  must  be  controlled  or  suppressed.  To  properly 
truncate  eq  (103),  the  left  side  will  have  the  top  left  12  x 
12  sub-matrix  extracted  from  each  45  x  45  matrix.  The  right 
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side  will  have  the  top  12  z  10  sab-matrix  extracted  from  the 
45  x  10  matrix  which  results  after  all  matrix  multiplication 
has  been  performed.  This  is  done  within  the  computer  program 
SHUTBM  (see  Appendix  C)  and  passed  to  the  ACOSS  program.  The 
smaller  matrix  equation  can  now  be  put  in  state-space  form. 

If  the  state  vector  is  defined  as 


x 


( 104) 


then  the  state  equation  can  be  written  as 


x  =  Ax  +  Bu 


( 105) 


where 


A  = 


B  = 


(106) 


The  matrix  Bf  will  depend  on  what  is  chosen  as  the  control 
vector  u(t).  For  this  analysis,  the  shuttle  can  be  torqued 
about  all  three  axes  as  can  the  antenna.  The  two  proof-mass 
actuators  can  each  produce  forces  in  the  X  and  Y  directions. 
Using  moments  and  forces  for  the  control  vector,  it  can  be 
written  as 


Fr2  Fp2  Fr3  Fp3  **lx  *4x  Mly  M4y  "li  «4r 
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(  107) 


F r 2  and  Fp2  are  the  forces  is  the  roll  and  pitch  axes 
directions  made  by  the  actuator  at  S=40. 

F  r  3  and  Fpj  are  the  forces  in  the  roll  and  pitch  axes 
directions  made  by  the  actnator  at  S=80. 

Mlx,  M^y>  and  M^z  are  the  moments  applied  to  the  shuttle. 

M  4  z »  M^y,  and  M^z  are  the  moments  applied  to  the  antenna. 

This  control  vector  must  be  factored  out  of  the  term  on 
the  right  side  of  eqn  (103)  before  it  is  truncated.  Vhat  is 
left  of  the  F  matrix  after  factoring  (and  before  premultipli¬ 
cation  by  the  E  matrix)  is  denoted  as  F^>  This  matrix  must 
be  modified  slightly  so  that  the  elements  of  the  control 
vector  will  be  roughly  of  the  same  magnitude.  Reference  [2] 
infers  that  the  actuator  forces  will  be  in  units  of  pounds 
(since  the  masses  being  moved  weigh  10  pounds  and  are  driven 
a  distance  of  only  one  foot),  or  perhaps  tens  of  pounds.  The 
moments,  however,  are  limited  by  the  NASA  paper  [2]  to  10,000 
ft-lbs.  If  these  moments  are  expressed  in  the  control  vector 
in  units  of  thousands  of  ft-lbs,  then  the  forces  and  moments 
will  be  of  roughly  the  same  order.  To  make  the  conversion, 
all  elements  of  the  Fj,  matrix  corresponding  to  a  moment  in 
the  control  vector  must  be  multiplied  by  1000.  This  has  been 
done  in  the  computer  program  MODEL  (see  Appendix  C). 

The  sensor  output  can  be  expressed  as 

y  -  C  x  (108) 


For  the  problem  at  hand,  the  output  vector  y  is  given  by 


y  -  j  e3  ur(sn2)  Op(sn2)  ur(sn3)  Up(sn3)  j  (109) 

The  ACOSS  program  takes  the  truncated  matrix  equation  and 
forms  the  state-space  equation.  It  also  must  be  given  the  C 
matrix  before  it  can  start  the  control  algorithm.  Once  it 
has  this  information  along  with  input  options,  it  can  begin 
forming  the  control  law,  a  brief  description  of  which  fol¬ 
lows. 

As  both  Janiszewski  [6]  and  Aldridge  [71  point  out,  the 
state  x  cannot  be  measured  directly.  It  can  only  be  measured 
through  the  output  y,  so  a  state  estimator  has  been  developed 
for  use  with  this  control  technique.  This  estimator  takes 
the  output  y  and  makes  a  best  estimate  of  the  state  x  which 
corresponds  to  y.  The  result  of  this  estimator  is  an  equa¬ 
tion  for  the  estimated  error: 

ec(t)  =  <AC  -  KcCc)  ec(t)  (110) 

f 

which  includes  the  estimator  gain  Kc*  which  was  found  through 
minimizing  the  quadratic  regulator  performance  index  [8:537]. 

The  same  process  is  applied  in  finding  the  control  gain  G 
to  be  used  for  feedback  purposes  in  the  equations: 

xc(t)  =  (Ac  +  B  CG ) x  c ( t )  +  BcGe(t)  (111) 

*s  =  +  BsGic(t)  +  BsGe(t) 

The  controller  and  observer  gain  matrices  E  and  G  are 
determined  such  that  the  performance  indices  J  and  are 


minimized  where 


c  -  J V« 

J0 


cxc  +  u  Rcu)  dt 


(112) 


J-  CO 

<«cT® 

0 


oec  +  vTrov)  dt 


The  matrices  Qc,  Qq,  Rc,  and  RQ  are  chosen  by  the  control 
designer.  If  Qc  and  Q0  are  chosen  positive  semi-definite  and 
Rc  and  RQ  positive  definite,  then  AC+BCG  and  Ac~KCc  are 
guaranteed  stable. 

The  closed  loop  system  model  can  now  be  formed,  as  Janis- 
zewski  [6]  by  defining  a  new  state  vector: 


z  (  t )  = 


-  T  1 

«.T<‘>  ! 


icT(t)  |  X,T(t) 


(113) 


which  incorporates  the  controlled  states,  suppressed  states, 
and  estimator  error.  The  closed  loop  system  model  can  then 
be  expressed  as: 


z  ( t )  = 


A  c  +  BcG 


i 


®cG 


r 


I  Ac  -  KC 


c  I  KCs 


-I- 


bsg 


I 


®sG 


I 


i(t)  (114) 


The  eigenvalues  of  the  above  matrix  will  show  the  stabi¬ 
lity  (or  instability)  of  the  system.  All  negative  eigen¬ 
values  (or  complex  conjugate  eigenvalues  with  a  negative  real 
part)  will  show  the  system  to  be  stable  within  the  limita¬ 
tions  of  the  model.  Any  positive  real  parts  of  eigenvalues 
will  show  the  system  to  be  unstable.  This  would  be  caused  by 
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a  coupling  effect  of  the  KC$  term  (called  observation  spill¬ 
over)  and/or  the  B ,G  terms  (called  control  spillover)  since 
the  optimal  regulator  theory  insures  that  AC+BCG  and  Ac-KCc 
are  stable  matrices. 

Once  the  ACOSS  program  has  built  eqn  (114),  it  runs  the 
eigenvalue  problem.  At  this  point  it  is  an  unsuppressed  run. 
The  program  will  then  run  the  suppression  algorithm,  which 
effectively  stabilizes  the  system  if  it  was  initially  un¬ 
stable. 

The  suppression  algorithm  can  be  implemented  in  one  of 
two  ways.  The  intent  is  to  drive  eqn  (114)  into  an  upper  or 
lower  diagonal  form  so  that  the  eigenvalues  of  the  matrix 
will  be  the  eigenvalues  of  the  terms  on  the  diagonal,  which 
all  have  a  negative  real  part.  This  can  be  done  by  either 
dr  iv ing 


B  8G  =  0 

or 

KCS  =  0 

Care  must  be  taken  to  keep  from  allowing  BCG  or  KCc  to  become 
zero,  lest  control  or  observation  be  completely  lost. 

The  method  used  by  ACOSS  is  to  find  a  transformation 
matrix  T  such  that  the  new  control  vector  U(t)  will  become 


u ( t )  =  Tv( t ) 


(115) 


This  is  done  through  a  technique  known  as  Singular  Value 
Decomposition  [91.  Aldridge  [71  gives  a  straightforward 
explanation  of  SVD  on  pages  42-48. 

After  having  found  the  T  matrix,  ACOSS  then  reruns  the 
eigenvalue  problem,  with  the  results  being  eigenvalues  with 
negative  real  parts.  The  system  is  thus  shown  to  be  stable 
within  the  bounds  of  the  mathematical  model. 

The  next  chapter  will  show  what  the  program  did  to  the 
shuttl e-beam- antenna  system.  It  is  meant  as  a  guide  to  show 
that  the  system  can  indeed  be  stabilized,  and  what  effect 
different  control  weightings  have  on  its  stability. 


t  f 


V .  Re  sal t  s 

The  intent  of  this  investigation  was  to  show  whether  or 
not  the  mathematical  model  developed  in  Chapter  III  coaid  be 
stabilized,  and  how  the  closed  loop  damping  could  be  improved 
over  the  open  loop  damping. 

No  time  response  was  calculated  for  this  system.  The 
measure  of  performance  was  chosen  to  be  the  closed  loop  modal 
damping  coefficients,  After  each  run,  this  is  compared 
to  each  open  loop  damping  coefficient,  (0£,  which  was  given 
(see  Table  I)  as  0.003.  A  target  value  for  improvement  was 
set  at  .03,  which  is  a  factor  of  10  above  the  open  loop 
coefficient. 

An  initial  run  of  the  ACOSS  program  was  made,  using  zero 
initial  conditions.  Initial  weighting  values  of  1  were  given 
to  the  weighting  matrix  F  of  the  natrix-Riccati  equation 
[8:541].  The  modes  selected  for  control  were  those  with  the 
six  lowest  frequencies.  These  were  the  three  rigid  body 
modes  and  the  three  lowest  roll,  pitch,  and  yaw  modes.  The 
suppressed  modes  were  those  four  with  the  next  lowest  fre¬ 
quencies,  and  the  residual  modes  were  the  remaining  two, 
which  had  the  highest  frequencies  of  those  modeled.  ACOSS 
ran  both  the  unsuppressed  and  the  suppressed  algorithms.  The 
results  are  shown  in  Table  IV  and  Table  V.  Table  IV  shows 


all  of  the  eigenvalues  of  the  modeled  system  to  have  negative 


real  parts.  Table  V  shows  that  the  .03  performance  index  was 
easily  surpassed  with  these  weightings.  There  is  very  little 
difference  between  the  nnsnppressed  and  suppressed  portions 
of  the  run,  mainly  because  the  system  was  stable  to  begin 
with,  and  no  suppression  was  needed. 

Another  run  was  made  with  the  weightings  on  the  control¬ 
led  flexible  modes  set  at  SO  to  determine  the  effect  on 
system  eigenvalues.  The  results,  shown  in  Table  VI  and  Table 
VII,  as  expected,  show  very  heavy  damping  on  the  flexible 
modes  and  a  very  slight  change  on  the  rigid  body  mode 
damping.  The  system  is  still  stable  without  suppression,  so 
the  suppression  portion  showed  little  change  as  before. 

Another  run  was  made  to  see  if  the  system  eigenvalues 
could  be  driven  unstable.  The  weighting  of  the  flexible  yaw 
mode  was  set  at  1000.  The  system  did  indeed  become  unstable, 
as  shown  by  the  overall  system  eigenvalues  in  Table  VIII. 
There  are  three  complex  conjugate  pairs  of  eigenvalues  with 
positive  real  parts  resulting  from  this  weighting.  The  sup¬ 
pression  algorithm  stabilized  the  system,  however,  as  shown 
in  Table  VIII  and  Table  IX.  The  eigenvalues  again  all  have 
negative  real  parts,  showing  the  system  to  be  stable.  This 
shows  that  the  ACOSS  program  works  on  this  model  as  it  was 


designed. 


Table  IV 


Overall  System  Eigenvalues 
Initial  Ran 


Controlled 

Mode  s:  1,2, 3, 4, 5, 6 

Suppressed 

Modes:  7,8,9,10 

Residual  Modes:  11,12 

Qc  =  Qo  =  M 

Rc  =  R0  =  M 

Before 

Sutpre s  s ion 

A_f _t r  Sunnr  e  s  s  ion 

-0.24024311 

+ 

80 .06607536  i 

-0 .24019691 

+ 

80 . 06608979  i 

-0  .24055247 

+ 

80 .18813210  i 

-0.24056484 

+ 

80 .18812917  i 

-0 .10488282 

+ 

34 . 81784403  i 

-0 .10445463 

+ 

34 . 81805332  i 

-0 .10461379 

+ 

34 .94614859  i 

-0 .10483  854 

+ 

34 .946 0227 4i 

-0 .03907  432 

+ 

13  .00359541  i 

-0 .03902331 

+ 

13  .00771146  i 

-0.03814158 

+ 

13  .13582349  i 

-0 .03942096 

+ 

13  .14026087  i 

-2 .39955811 

+ 

5 .086  87  8294  i 

-2 .40027385 

+ 

5 .076907440  i 

-0.16797683 

+ 

5.594573741 i 

-0.01679353 

+ 

5  .5978098091 

-0 .76047227 

+ 

1  .785963613  i 

-0 .75564842 

+ 

1  .762731235  i 

-0.70649699 

+ 

1 .803664617  i 

-0 .70786680 

+ 

1  .776586  936  i 

-0 .157207  56 

+ 

1  .  845696247  i 

-0.00556498 

+ 

1  .  854871652  i 

-0.15640709 

+ 

1  .836151481  i 

-0.00553401 

+ 

1 .844585698  i 

-0 . 86594036 

+ 

0 .499810091  i 

-0 . 86602637 

+ 

0 .4999948081 

-0.86622385 

+ 

0 .500038680  i 

-0.86603015 

+ 

0 .4  99998044 i 

-0 .02637183 

+ 

0 .026597291  i 

-0 .02661848 

+ 

0 .0266011  83  i 

-0.00990540 

+ 

0.009912217 i 

-0.00991154 

+ 

0.009910617 i 

-0 .00994017 

+ 

0 .010025435  i 

-0 .01002882 

+ 

0 .010027  879 i 

-0.86602540 

+ 

0 .500000000  i 

-0 .86602540 

+ 

0  .499999999  i 

Table  V 

Eigenvalues  and  Closed  Loop  Damping  Coefficients 


Initial  Rnn 


Controlled  Modes:  1,2, 3, 4, 5, 6 
Suppressed  Modes:  7.8,9.10 
Residual  Modes:  11,12 

Qc  =  fto  =  M 

Kc  -  Ro  -  M 

Controlled  Modes: 

Before  Sunn r ession 


After  Sui 


ression 


-0.00991 

+ 

0.00991  i 

5 

a 

.707 

-0.00991 

+ 

0.00991 i 

5 

.707 

-0 .01003 

+ 

0.01003  i 

5 

= 

.707 

-0.01003 

+ 

0  .01003  i 

5 

= 

.707 

-0 .02662 

+ 

0 . 0  26  6  0  i 

5 

= 

.707 

-0.02662 

+ 

0 .02660  i 

5 

= 

.707 

-2.40027 

+ 

5  .07691  i 

5 

= 

.427 

-2.00427 

+ 

5  .07691  i 

5 

= 

.427 

-0.707  87 

+ 

1 .77659  i 

5 

= 

.370 

-0 .707  87 

+ 

1 .77659  i 

t 

= 

.370 

-0.75565 

+ 

1 .76273  i 

5 

= 

.394 

-0.75565 

+ 

1 .76273  i 

5 

= 

.394 

Suppressed  Modes: 


Before  Sunression 


-0  .03902 

+ 

13  .0077  i 

?  = 

.003 

-0 .03902 

+ 

13  .0077  i 

5  = 

.003 

-0  .03942 

+ 

13 .1403  i 

5  - 

.003 

-0.03942 

+ 

13 .1403  i 

5  = 

.003 

-0.10445 

+ 

34 .8181  i 

5  = 

.003 

-0.10445 

+ 

34.8181 i 

5  = 

.003 

-0.10484 

+ 

34 .9460i 

5  = 

.003 

-0.10484 

+ 

3  4 . 946  0 i 

5  - 

.003 

Residual  Modes: 

Before  Suppression 

-0.24019  +  80  .0661  i  $  =  .003 
-0.24056  +  80 . 1 8 81  i  5  =  .003 


-0.24019  ±  80 . 0 6 6 1  i  £ 
-0.24056  +  80  .1881  i  5 


X 


Controlled  Modes:  1,2, 3, 4,5, 6 


Suppressed  Modes:  7,8,9,10 
Residual  Modes:  11,12 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Rc  =  Rc  =  W 


0  0  0  0 
10  0  0 
0  10  0 
0  0  50  0 

0  0  0  50 

0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 


0  0 
0  0 
0  0 
0  0 
0  0 
50  0 
0  1 
0  0 
0  0 
0  0 
0  0 
0  0 


0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 


0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
50  0  0 

0  50  0 

0  0  50 


Before  Suppression 


After  Suppression 


-0.24000054  ± 
-0  .23968480  + 
-32.9177464  + 
-0.09589029  + 
-0  .09286895  + 
-0.07403583  + 
-0.06884054  + 
-8.07978974  + 
-7.47846553  + 
-0.27405494  ± 
-0.55446360  + 
-0.53340846  ± 
-1.48750669  + 
-0.66519955  + 
-0.50217705  + 
-0.02355995  + 
-0.00994269  + 
-0.00877420  + 
-1 .94580388  ± 
-1.94524442  ± 
-1.94564483  + 


80 .06585026  i 
80  .18829557  i 
0  i 

34.81961022  i 
34 .951973  91  i 
12.88756352  i 
13  .01455579  i 
0  i 
0  i 

5 .586351641  i 
1  .  858837822  i 
1.819136002  i 
0  i 
0  i 
0  i 

0.025848884  i 
0 .009989706  i 
0.009894585  i 
1 . 813202193  i 
1 .812357  53  4i 
1  .  812604131  i 


-0.24018449 
-0.24056131 
-32 . 9345667 
-0.10445463 
-0.10483854 
-0  .03902331 
-0 .03942096 
-8.43563283 
-7 . 96961656 
-0.01679357 
-1 .11786810 
-0  .00556567 
-0 .00553445 
-0 .39513958 
-0 . 43372921 
-0 .02705426 
-0 .00993458 
-0 .01005245 
-1 .94563452 
-1 .94563911 
-1.94564485 


+  80.06609121i 
+  80 .1881296  8i 
+  0  i 

+  34.818053311 
+  34 .94602274i 
+  13 .00771146  i 
±  13 .140  260  87  i 
+  Oi 
+  0  i 

+  5 .597  809809  i 
+  0  i 

±  1 .854871646  i 
+  1 . 844585695  i 
+  Oi 
+  0  i 

+  0 .02616097  8i 
+  0 .009887539  i 
±  0 .010004210  i 
+  1  .  812601329  i 
+  1 .812608398  i 
+  1  .  812604177  i 


Eigenvalues  and  Closed  Loop  Damping  Coefficients 


Controlled  Modes:  1,2, 3,4 ,5,6 
Suppressed  Modes:  7.8,9,10 
Residual  Modes:  11,12 
Weightings:  See  Table  VI 

Controlled  Modes: 


ore  Sunoressioi 

-0.00993  +  0  .00989 i  $ 
-0  .01005  +  0  .01000 i  5 
-0  .02705  +  0  .02616 i  $ 
-0.39514  +  Oi  ? 

-0.43373  +  Oi  $ 

-1  .11787  +  0  i  5 

-7  .96596  +  0  i  5 

-8.43561  +  Oi  5 

-32.9346  +  Oi  ? 


.709  -0.00993  +  0.00989i 

.709  -0.01005  ±  O.OlOOOi 

.719  -0  .02705  +  0  .02616  i 

1  .0  -0.39514  +  Oi 

1  .0  -0.43373  +  0  i 

1  .0  -1  .11787  +  Oi 

1 .0  -7.96596  +  0  i 

1  .0  -8.43561  +  Oi 

1.0  -32.9346  +  Oi 


.709 
.709 
.719 
1  .0 
1  .0 
1  .0 
1  .0 
1  .0 
1  .0 


Suppressed  Modes: 

Before  Sunpressioi 

-0  .03902  +  13  .0077  i  $ 
-0.03942  +  13.1403 i  $ 
-0.10445  +  34.8181 i  $ 
-0.10483  +  34 ,9460i  $ 

Residual  Modes: 


003 

-0  .03902 

+ 

13  .0077  i 

5  = 

.003 

003 

-0.03942 

+ 

13 .1403  i 

5  = 

.003 

003 

-0.10445 

+ 

34.8181 i 

5  = 

.003 

003 

-0.10483 

+ 

34 .9460i 

5  = 

.003 

Before  Sunnre 

-0.24019  ±  80  .0661  i 
-0.24057  +  80 . 1 881  i 


■0.24019  ±  80 .0661  i 
-0.24057  +  80 . 1  881  i 


Controlled  Modes:  1,2, 3, 4, 5, 6 


Suppressed  Modes:  7,8,9,10 


Residual  Modes:  11,12 


-150.26995172  +  Oi 
-0  .2411541224  +  80.0589083  1 
-0.2059402063  +  80.2074037i 
-38.280270629  +  Oi 
-0.0865376915  +  34.80520091 
+0.1107479215  +  35.2931502i 
+0.0272984940  +  12.8537225i 
+0.5538430430  +  11  .0878564i 
-7.3  933793757  +  Oi 
-0.9721057460  +  5.50363776i 
-2.8091057926  +  1.325751411 
-0.5164799305  ±  1.738385781 
-1  .3225958070  +  0  i 
-0.4337955132  +  Oi 
-0.2052694774  +  0  i 
-0.0224186944  +  0.02557039i 
-0.0091808952  +  0.00986227  i 
-0.0070312672  +  0.00942698i 
-1.9451667499  +  1.81280739i 
-1.9457737116  +  1.81253073i 
-1.9456447513  +  1.81260413i 


-150.300112  +  Oi 
-0.24018206  +  80 .0660909 i 
-0.24056305  +  80.1881308i 
-39.0747709  +  Oi 
-0.10445463  +  34.8180533i 
-0.10483854  +  34.9460227  1 
-0.03902331  +  13 . 0077115 i 
-0.03942096  +  13.1402609i 
-8.19451804  +  Oi 
-0.01679430  +  5.5978098H 
-0.00558402  +  1.85487153i 
-0.00553445  +  1.84458569i 
-0.45336639  +  Oi 
-0.21174024  +  Oi 
-0.09184072  +  Oi 
-0.02703085  +  0.026180731 
-0.01003086  +  0 . 010023  82  i 
-0.00992552  +  0.00989796  i 
-1.94563263  +  1.81260802i 
-1.94564374  +  1.812604811 
-1.94564486  +  1.812604181 
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Eigenvalues  and  Damping  Coefficients 


Controlled  Modes:  1,2, 3, 4. 5. 6 
Suppressed  Modes:  7,8,9,10 
Residual  Modes:  11,12 


Weightings 

:  See  Table 

VII 

Before  Sunnression 

After 

Sunnre  s  s ion 

-0  .00992 

+ 

0 .00998i 

5 

= 

.708 

-0.00992 

+ 

0 . 00  99  8  i 

5 

= 

.708 

-0.01003 

+ 

0 .01002i 

5 

= 

.707 

-0  .01003 

+ 

0  .01002 i 

5 

cs 

.707 

-0.02731 

+ 

0 .02618i 

i 

= 

.718 

-0.02731 

+ 

0  . 026 1 8 i 

5 

z= 

.718 

-0  .09184 

+ 

0  i 

5 

= 

1  .0 

-0.09184 

+ 

0  i 

i 

= 

1  .0 

-0.21174 

+ 

0  i 

5 

a 

1  .0 

-0.21174 

+ 

0  i 

= 

1  .0 

-0  .45336 

+ 

0  i 

5 

= 

1.0 

-0  .45336 

+ 

0  i 

= 

1  .0 

-8.19445 

+ 

Oi 

$ 

a 

1  .0 

-8.19445 

+ 

0  i 

5 

= 

1.0 

-39  .0747 

+ 

Oi 

$ 

= 

1.0 

-39.0747 

+ 

Oi 

= 

1  .0 

-150.300 

+ 

0  i 

$ 

= 

1  .0 

-150 .300 

+ 

0  i 

5 

= 

1.0 

Suppre  s  sed 

Modes: 

Before  Sunnre 

s s  ion 

After 

Sunnr e  s  s ion 

-0  .03902 

+ 

13  .0077  i 

5 

— 

.003 

-0  .03902 

+ 

13  .0077  i 

.003 

-0  .03942 

+ 

13  .1403  i 

5 

= 

.003 

-0.03942 

± 

13.1403  i 

5 

= 

.003 

-0 .10445 

+ 

34 .8181 i 

= 

.003 

-0 .10445 

+ 

34.8181 i 

5 

= 

.003 

-0.10483 

+ 

34 .9460i 

* 

.003 

-0.10483 

+ 

34 .9460  i 

S 

= 

.003 

Residual 

Modes: 

Before  Sunnression 

After 

Sunnr e  s  s 

ion 

-0.24019 

+ 

80  .0661  i 

5 

= 

.003 

-0.24019 

+ 

80  .0661  i 

5 

__ 

.003 

-0.24056 

+ 

80  .1881  i 

5 

.003 

-0 .24056 

+ 

80 .1881  i 

5 

= 

.003 

V .  Conclusions 

A  mathematical  model  was  developed  for  a  shuttle-beam- 
ante  nna  system.  The  system  was  discretized  using  an  assumed 
modes  approximation.  The  equations  of  motion  were  developed 
from  scratch  and  linearized  by  ignoring  higher-order  terms 
and  assuming  small  beam  deflections.  The  coupling  between 
the  rigid-body  and  flexible  motions  was  taken  into  account  in 
the  equations  of  motion.  Fourteen  modes  were  assumed  in  each 
of  the  three  directions  of  motion  (roll  and  pitch  bending  and 
yaw  torsion).  The  equations  of  motion  were  put  in  matrix 
form,  and  the  matrices  diagonalized. 

Linear  optimal  regulator  theory  was  applied  to  examine 
the  stability  of  the  mathematical  model.  A  target  value  of 
.03  for  the  closed  loop  modal  damping  coefficient  was  used  as 
a  measure  of  improvement  between  the  open  loop  and  closed 
loop  system.  The  target  value  was  surpassed  on  the  first 
try,  using  equal  weightings  on  all  of  the  modeled  modes.  The 
weightings  were  then  modified  to  examine  the  eigenvalues  of 
the  system  for  open  loop  instability.  This  was  found  with  a 
very  large  weighting  of  the  flexible  yaw  mode.  The  closed 


loop  suppression,  as  expected,  r e- s t ab i 1 iz ed  the  system. 


VII. 


Re commend at  ions 

The  major  emphasis  of  this  analysis  was  on  the  mathemati¬ 
cal  modeling  of  the  shut t 1 e-b e am- ant e nna  system.  Unfor¬ 
tunately.  only  a  cursory  investigation  of  the  control  of  the 
system  could  be  accomplished.  A  complete  study  should  be 
made  of  the  time  response  of  the  system,  using  varying  ini¬ 
tial  conditions  to  measure  how  quickly  the  system  can  be 
controlled. 

As  mentioned  in  Appendix  A.  the  choice  for  the  locations 
of  the  two  actuators  was  based  solely  upon  the  mode  shapes  of 
a  cantilever  beam  with  no  mass  attached  to  its  free  end.  An 
analysis  could  be  made  to  determine  the  best  locations  for 
the  actuators,  incorporating  the  mass  of  the  antenna  into  the 
mode  shapes  of  the  beam. 

The  NASA  Challenge  [2]  also  gave  information  (which  needs 
to  be  carefully  corrected)  on  1 i n e- o f- s i g h t  calculations. 
This  would  be  a  better  measure  of  system  performance  than  the 
closed  loop  damping  coefficient  and  should  be  the  subject  of 
further  study.  The  paper  also  sets  up  a  maneuvering  problem, 
incorporating  attitude  changes  and  a  slew  maneuver  of  the 
antenna.  A  complete  study  can  be  made  on  this  subject  alone. 

The  major  building  blocks  of  a  thorough  investigation  of 
the  control  of  the  system  have  now  been  created.  The  compu¬ 
ter  programs  have  been  built  with  the  flexibility  to  alter 


t. 


< / 


3. 


5 

6. 
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Cantilever  Mode  Shapes 

The  NASA  paper  [2]  gives  graphs  for  the  roll,  pitch,  and 
yaw  mode  shapes  of  the  beam  without  reference  to  the  type 
model  used.  An  investigation  into  a  fixed-free  cantilever 
beam  and  a  free-free  beam  showed  that  the  mode  shapes  of  a 
fixed-free  beam  closely  resemble  those  in  the  NASA  paper. 
The  associated  frequencies  do  not  match  exactly,  but  no  clue 
was  given  as  to  the  number  of  modes  assumed  by  the  paper,  and 
therefore  the  accuracy  of  the  frequencies  cannot  be  measured. 
The  next  few  pages  show  the  first  seven  calculated  mode 
shapes  for  a  fixed-free  cantilever  beam  using  14  modes  for 
the  roll  (and  pitch)  motion,  14  modes  for  the  yaw  motion. 
Since  the  beam  bending  is  identical  for  both  roll  and  pitch, 
the  roll  graphs  also  accurately  depict  the  pitch  modes.  A 
comparison  with  the  graphs  in  the  NASA  paper  shows  the  close 
re  sembl ance . 

It  was  necessary  to  choose  the  positions  along  the  beam 
for  the  locations  of  the  two  proof-mass  actuators.  This  was 
done  by  inspection  of  the  roll  (and  therefore  pitch)  mode 
shapes  using  the  14-mode  approximation  graph.  The  higher  the 
amplitude  of  any  mode  shape  at  a  point  along  the  beam,  the 
more  that  particular  mode  will  be  influenced  by  an  actuator 
located  at  that  point.  A  best  fit  was  thus  made  to  the  seven 
curves  inspected.  The  40-foot  point  showed  adequate  ampli- 


tude  for  all  seven  modes  and  was  chosen  as  the  position  of 
actuator  number  one.  The  80-foot  point  showed  adequate  am¬ 
plitude  for  all  but  mode  number  5.  Since  the  frequency  for 
mode  5  is  fairly  high  and  adequate  amplitude  was  shown  at  the 
40-foot  point,  it  was  decided  that  the  80-foot  point  would  be 
adequate  for  the  position  of  actuator  number  two.  It  should 
be  kept  in  mind  that  these  are  the  mode  shapes  for  a  canti¬ 
lever  beam  without  any  mass  at  its  free  end.  Since  the 
s h u t t 1 e- b e a m- a n t e n n a  system  does  indeed  have  such  a  condi¬ 
tion,  the  best  positions  for  the  actuator  locations  could 
very  well  be  different  from  those  used  here.  This  would 
represent  a  completely  different  study,  and  the  positions 
chosen  here  will  serve  the  purposes  of  this  investigation. 
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Fig.  8.  Yaw  Mode  Shapes 


D.e  r_i  vji_t_ij5  n  o_f  Coefficients  for  Roll  and  P  j._t_ch  Fane  t  ions 

This  section  shows  the  development  used  to  solve  for  the 
coefficients  of  the  equation 

Ri  =  A^sinaL  +  B  •  cosaL  +  C^sinhaL  +  D- coshaL 

with  the  constraint 

fL 

1  p  AR  ^R ^  ds  =  1 

J  0 

The  i  subscript  will  be  dropped  for  convenience,  and  the  p  A 
will  be  rewritten  as  pa  to  avoid  confusion.  Several  other 
constraints  also  apply  to  this  problem: 

A  =  -C 

B  =  -D 

cosaL*coshaL  =  -1 

B  =  -A(sinSL  +  sinhaL)  =  - Q A 
(cosaL  +  coshaL) 

Thus  the  function  becomes 


R  =  AsinaL  -  QAcosaL  -  AsinhaL  +  QAcoshaL 


Substituting  this  function  into  the  constraint  equation 

f  L 

1  =  pa  |  [Asinas  -  QAcosas  -  Asinhas  +  QAcoshas]^  ds 


QAco  s  a  s 


Carrying  oat  the  multiplication  results  in 
1/pa  =  A^j‘^[sin^os  -  2Qsinas  -  2sinassinhos 


+  2Qsi nas co shas  +  Q2cos2as  +  2Qco sas s i nhas 


-2Q^ co sas co shas  +  sinh2as  -  2Qs i nhas co shas 
+  Q^cosh^as]  ds 

Integrating  each  term  separately  yields 

1/pa  =  -  sin2as  1  L  -  2Q  s i n2a s  ^ 

i  I  2  4a  0  2a  0 


-  2  coshassinas  ~  sinhascosas 

L  2a  Jo 

+  2Q  sinhassinas  -  co shas co sas I ^ 

2a  0 


+  Q2  s  +  sin2a  s  ^ 


4a  0 


+  2Q  coshascosas  +  sinhassinas 

2a  0 

2  1  f 

-  2Q  sinhascosas  +  coshassinas 

2  a  0 


sinhascoshas  -  s 

L  -  2Q 

s i nh2a  s 

2a  2 

0 

2a 

s i nhas  co  shas 


Plugging  in  the  limits  and  rearranging  terms  gives: 


1/pa  =  A2 


-  1  |  sin2aL  -  Qsin^aL 
a 


4a 


(sirLi) 


coshaLsi naL 


sinhaLcosaL 


+  2QsinhaL  sinaL  +  |  Q2  +  1  \  sinhaLcoshaL 
a  \  2a  / 


+  Q2L  -  Qsinh2aL 


To  simplify  this  further,  a  very  close  look  must  be  taken  at 


Q,  Q2.  Q2  -  1,  and  Q2  +  1 : 


Q2  = 


(sinaL  +  sinhaLy 
cosaL  +  co shaL  I 


=  sin2aL  +  2sinaLsinhaL  +  sinh2aL 


cos2aL  +  2cosaLcoshaL  +  cosh^aL 


but 

cos2aL  +  2cosaLcoshaL  +  cosh2aL 
=  cos2aL  -2  +  cosh2aL 
=  (cos2aL  -  1)  +  ( co  sh2aL  -  1) 

=-sin2aL  +  sinh2aL 


so 


substituting  : 


<5 


Q^-l  =  sinhaL  +  2 s i naL s i nhaL  +  sinh^aL  -  sinhaL  +  sinh^aL 


-sinzaL  +  sinhzaL 


=  2sin^aL  +  2 s i naL s i nhaL 


-sinhaL  +  sinhzaL 


in  the  same  manner: 


+  i  =  2sinaLsinhaL  +  2sinh^aL 


-sinhaL  +  sinhzaL 


Substituting  in  these  expressions  for  Q,  Q^,  Q^-l,  and 


1  /  p  a  =  A 


=  *2 


si naL/  si naL  +  sinhsL  \  sin2aL 


2a  V-sin^aL  +  sinh^sL 


(sinaL  +  s inhaL ) ( s i n^qL ) 
a(cosaL  +  coshaL) 


-2  ( sinaLsinhaLcoshaL) ( sinaL  +  sinhaL) 


(-sinzaL  +  sinh^aL) 


-  2$in^aLcosaL( sinaL  +  sinhaL) 


a(-sin^aL  +  sinh^aL) 


+  2  s  i  naL  s  i  nhaL  (  s  i  naL  +  sinhaL) 


a(cosaL  +  coshaL) 


+  s i nh^ aL co  shaL (si naL  +  sinhaL) 
a(-sin^aL  +  sinh^aL) 


+  L(sinzaL  +  2 s i naL s i nhaL  +  sinh^aL) 


now,  let 


so  : 

1/  p  a  = 


+ 


+ 


Comb ining 
1/pa  = 


+ 


-sinhaL  +  sinh^aL 


s i nh^aL ( s i naL  +  sinhaL) 
a ( co  saL  +  coshaL) 


si naL  +  sinhaL  =  z 


-sinhaL  +  sinh^aL 


7.  7 

isin^aLcosaL  -  zsin4aL( cosaL  +  coshaL) 
"a  a 


2xs i naL s i nhaL co shaL  -  2 x s i naL co s aL s i nhaL 
a  a 


2x(cosaL  +  co shaL ) s inhaL si naL  +  x s i nh^aL co shaL 
a  a 


L/sinaL  +  sinhaL^ 


cosaL  +  coshaL 


x_(cosaL  +  coshaL)  sinh^aL 
a 


and  cancelling  appropriate  terms: 


-_xs  i  n^aLco  shaL  -  xs  inh^aLco  saL 
a  a 


L /sinaL  +  sinhaL! ^ 
IcosaL  +  coshaL/ 


Combining  the  x/a  terms  and  r e- sub s t i t n t i ng  for  x: 


So,  the  remaining  term  is  the  only  non-zero  term,  resulting 
in  : 


1 / p a  =  L  / sinaL  +  sinhaL\ ^ 

l  co  saL  +  co  shaL / 


And  solving  for  B: 


n  o  o  n 


PROGRAM  MODEL 


C 

C  . . . 

c 

C  THIS  PROGRAM  TAKES  THE  ANGLE  DATA  FROM  ' FDATA1 '  AND  COMPUTES 

C  THE  A- ,  B- ,  AND  C-  COEFFICIENTS  FOR  THE  ROLL,  PITCH,  AND  YAW 
C  EQUATIONS.  IT  ALSO  COMPUTES  THE  OMEG A-  SQUARES  FOR  THESE  EQUA- 
C  TIONS.  IT  THEN  EVALUATES  THESE  EQUATIONS  AT  BOTH  ENDS  OF  THE 
C  FIXED-FREE  BEAM.  IT  INCORPORATES  SOME  CORRECTION  FACTORS  DUE 
C  TO  THE  MOTION  OF  THE  PROOF-MASS  ACTUATORS  (LOCATED  AT  POSITIONS 
C  SN1  AND  SN2 ) .  THE  PROGRAM  PUTS  ALL  OF  THIS  INFORMATION  IN  MATRIX 

C  FORM,  CREATING  TWO  MATRICES — A  MASS  MATRIX  'M'  (STORED  IN  DATA 

C  FILE  ' FMASTM ’ ) ,  AND  A  STIFFNESS  MATRIX  'O'  (MADE  UP  OF  THE  OMEGA- 

C  SQUARES  AND  STORED  IN  DATA  FILE  ' FMASTK ' ) .  THE  REST  OF  THE  PRO- 

C  GRAM  SOLVES  THE  EIGENVALUE  PROBLEM  FOR  THESE  TWO  MATRICES  FOR 
C  COMPARISON  WITH  THE  NASA  PAPER. 

C 

c 

DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
C 

DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
C 

REAL  AA( 42 .42) , BB( 42 ,42) , COMP1 , COMP2 , VREAL( 42 ,42) , HERTZ ( 42) 
REAL  WK( 1888) , BBETA( 42 ) . FBMATX( 45 , 10) , PMATX(7 , 45) 

C 

COMPLEX  ALFA (42) ,VEC(42,42) ,VECT(42,42) ,EIGENV(42) .OMEGA (42) 
COMPLEX  TEMP (42) 

INTEGER  I.  J,  K.  L.  Q,  U,  V.  W.  X,  Y,  Z 
INTEGER  IA,  IB.  N,  IJOB,  IZ.  IER 
C 

PARAMETER  ( A1  =  90 5443. 00, A4  =  4969. 0.B1 =6789100. 0.B4  =  4 96  9.0, 

1  Cl  =  70 866  01  .0, C4  =  9  93  8 . 0 , D1 =-1453  93. 0, Ml  =  6366. 46, M4  =  12 .42 , 

1  M2 =. 3108, M3=. 3108, SN1=40 .0, SN2=80 .0) 

PARAMETER  ( V= 14 . W= 1 5 , X= 2 8 , Y= 2 9 , Z=4 2 ) 


A(14)  ,B(14)  ,C(14)  , ALPH  AL ( 14)  .ALPHA (14) 

BETAL( 1 4 ) , BETA ( 14) ,OMESQR(14) ,OMESQP(14) 

PH  I (14)  , PH IPR ( 14 )  , PHIZ (14)  ,PHIPRZ(14)  ,OMESQY(14) 
THETA ( 1 4 ) , THEPR (14) , THETAZ (14) , THEPRZ(14) 

PSI(14) , PS IPR (14) ,PSIZ(14) ,PSIPRZ(14) 

PHICR1 (14) , PH I CR2 (14) , THECR1 (14) . THECR2(14) 
PHICOR( 14 ,14) , THECOR (14,14) ,PMAT(7,45) 

FR ( 14 , 1 4 ) ,GR(14,14) ,FP(14,14) ,GP(14,14) .ONE 
FY( 14. 14) ,GY( 14,14) ,M( 42, 42) ,0(42,42) 

A1  ,  A4  .  B1  ,  B4  ,  Cl .  C4  ,  D1 ,  Ml ,  M4  ,  Y4  ,  Z4  ,  M2  ,  M3  ,  SN1  ,  SN2 
FB  MAT (45 , 10> 


INITIALIZATIONS 


Y4=- 1 8 . 7  5  *M4 
Z4=-32.5*M4 
ONE=- 1 .0 


n  n 


cccccccccccccccccccccccccccccccccccccccc 

C  READ  IN  THE  ANGLE  INFORMATION: 

cccccccccccccccccccccccccccccccccccccccc 

c 

OPEN ( UN IT=  7 , FILE= ' FDATA1 ' . ACCESS= ' SEQUENTIAL ' , STATU S= ' OLD' ) 

REW IND ( UN IT=7 ) 

10  FORMAT (E31 .24) 

DO  20  1=1 , 14 
READ (7.10) ALPH  AL ( I ) 

B( I ) =- 1/ ( SQRT( .09556*130)) 

A(I)=(-B(I))*(COS( ALPH AL (I))+COSH(ALPHAL(I)) ) / ( + SIN ( ALPH AL ( I ) ) 
1  +S INH( ALPHAL( I ) ) ) 

OMESQR(  I )  =  (  4E7/  (  .09556*130**4)  )  •  ( ALPH  AL  (  I )  *  *  4 ) 
OMESQP(I)=(4E7/( .09556*130**4) ) * ( ALPHAL( I ) **4) 

20  CONTINUE 
C 


C 


IF  (I.EQ.lll)  THEN 

PRINT*, 'YOUR  COEFFICIENTS  ARE:' 
PRINT*,'  ' 

DO  40  1=1,14 

WRITE (*,10) (A( I) ) 


40  CONTINUE 
C 

ENDIF 


DO  50  1=1,14 

BETAL( I)=( 2. 0*1-1) • (DACOS(ONE) ) /2 .0 
BETA ( I ) =BETAL ( I ) / 1 30 
50  CONTINUE 

DO  60  1=1,14 

C ( I ) =SQRT( 2/ ( .9089*130) ) 

OMESQY(I) =(4E7/ ( .9089*130**2) ) • ( BETAL ( I ) * • 2 ) 
60  CONTINUE 


IF(I.EQ.lll)  THEN 
C 

PRINT*,  'THE  BETA-LS  ARE:' 

DO  70  1  =  1 ,14 

WRITE ( * , 1 0 ) BETAL { I ) 

70  CONTINUE 

PRINT*, 'THE  BETAS  ARE:' 

DO  80  1=1,14 

WRITE(*, 10) BETA(I ) 

80  CONTINUE 

C 

PRINT*, 'AND  ALL  THE  CnS  ARE:' 
WRITE (*,10)C(1) 

C 


ENDIF 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  CALCULATE  THE  PHI,  PHI- PRIME,  AND  PHI  CORRECTION  FACTORS: 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c 

DO  90  1=1,14 

ALPHA(  I)  =  ALPHAL(  I)  /130 
90  CONTINUE 

DO  100  1=1,14 

PH  I  ( I ) = A ( I ) *  S IN ( ALPH AL ( I ) ) +B ( I ) *COS ( ALPH AL ( I ) ) 

1  -A(I)*SINH( ALPH AL ( I ) ) - B ( I ) *  COSH ( ALPH AL ( I ) ) 

C 

PH ICR1(I)=A(I)*SIN( ALPHA( I ) *SN1 )  +  B( I ) *COS( ALPHA( I ) *SN1 ) 

1  -  A ( I ) *S INH ( ALPHA (I)*SN1)  -  B ( I ) * COSH ( ALPHA ( I ) * SN1 ) 

C 

PHICR2 (I)=A(I)*SIN(ALPHA(I)*SN2)  +  B( I ) *COS ( ALPHA ( I ) *SN2 ) 

1  -  A(I) *SINH(ALPHA( I) *SN2)  -  B ( I ) *COSH( ALPHA ( I ) * SN2 ) 

C 

100  CONTINUE 

DO  110  1=1,14 

PHIPR(I)=( ALPH A ( I ) )*(A(I)*COS(ALPHAL(I) ) -B( I ) *S IN ( ALPH AL ( I ) ) 
1  —  A ( I ) *  COSH( ALPH AL ( I))-B(I)*SINH( ALPH AL ( I ) )  ) 

110  CONTINUE 

DO  120  1=1,14 
PH  12  (  I )  =0 
120  CONTINUE 

DO  130  1=1.14 
PH IPRZ ( I ) =0 
130  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  CALCULATE  THE  THETA,  THETA-PRIME,  PSI,  PSI-PRIME,  AND  THETA- 
C  CORRECTION  FACTORS: 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

DO  140  1=1,14 

THETA (I)=A(I)*SIN( ALPH AL (I))+B(I)*COS( ALPH AL ( I ) ) 

1  -A(I)*SINH( ALPH AL ( I ) ) - B ( I ) *  COSH ( ALPH AL ( I ) ) 

C 

THECR1 (I)=A(I)*SIN( ALPHA( I ) *SN1 )  +  B ( I ) *COS ( ALPH A ( I ) • SN1 ) 

1  -A(I)*SINH(ALPHA(I)*SN1)  -  B ( I ) • COSH ( ALPHA ( I ) * SNl ) 

C 

THECR2 (I)=A(I)*SIN(ALPHA(I)*SN2)  +  B ( I ) *COS ( ALPH A ( I ) *SN2 ) 

1  -A( I) *SINH(ALPHA( I) *  SN2 )  -  B ( I ) * COSH ( ALPHA ( I ) * SN2 ) 

C 

140  CONTINUE 

DO  150  1=1 ,14 

THEPR ( I )  =  ( ALPH  A(I) )*(A(I)*COS(ALPHAL( I) )-B( I ) • S IN ( ALPH  AL ( I ) ) 
1  -A(I)*COSH(ALPHAL(I) )-B(I)*SINH( ALPHAL( I ) ) ) 

150  CONTINUE 

DO  160  1  =  1  ,14 
THETAZ ( I ) =0 


CONTINUE 
DO  170  1=1,14 
THEPRZ ( I ) =0 
CONTINUE 
DO  180  1=1,14 

PS  I ( I ) =C  < I ) *  S IN ( BETAL ( I ) ) 

CONTINUE 
DO  190  1=1,14 

PS  IPR  ( I )  =  ( BETA (I))*COS(BETAL(I) )*C(I) 

CONTINUE 
DO  200  1=1,14 
PSIZ(I)=0 
CONTINUE 
DO  210  1=1,14 

PS IPRZ ( I ) =C ( I ) *B  ETA ( I ) 

CONTINUE 

CHECK  FOR  PROPER  CALCULATIONS 

IF  (I.EQ.lll)  THEN 

PRINT*,  'THE  PHI  FUNCTIONS  EVALUATED  AT  ZERO  ARE 
DO  220  1=1,14 

WRITE( *. 10) PHIZ( I) 

CONTINUE 
PRINT*.  ' 

BBI*ar§,  lOElfHI  FUNCTIONS  AT  130  ARE:' 

WRITE (*,10) PH 1(1) 

CONTINUE 
PRINT*, '  ' 

PRINT*. 'THE  PHI-PRIME  FUNCTIONS  AT  ZERO  ARE:' 

DO  240  1=1,14 

WRITE(*,10) PH IPRZ ( I ) 

CONTINUE 
PRINT*, '  ' 

PRINT*, 'THE  PHI-PRIME  FUNCTIONS  AT  130  ARE:' 

DO  250  1  =  1 .14 

WRITE(*,10)PHIPR(I) 

CONTINUE 
PRINT*, '  ' 


PRINT*,  'THE  THETA  FUNCTIONS  EVALUATED  AT  ZERO  ARE 
DO  260  1=1,14 

WRITE( *, 10) THETAZ ( I ) 

CONTINUE 
PRINT*. ' 

PRINT*, 'THE  THETA  FUNCTIONS  AT  130  ARE:’ 

DO  270  1=1,14 

WRITE ( * , 1 0 ) THETA ( I ) 

CONTINUE 
PRINT*.  '  ' 


290 

C 

300 


PRINT*, 'THE  THETA-PRIME  FUNCTIONS  AT  ZERO  ARE:' 

DO  280  1=1,14 

WRITE ( * , 10) THEPRZ( I) 

280  CONTINUE 

PRINT*. '  ' 

PRINT*, 'THE  THETA  PRIME  FUNCTIONS  AT  130  ARE:' 

DO  2  90  1  =  1 ,14 

WRITE(*, 10) TH  EPR ( I ) 

CONTINUE 
PRINT*, '  ' 

PRINT*, 'THE  ASSOCIATED  PSI  FUNCTIONS  AT  ZERO  ARE:' 
DO  300  1=1 ,14 

WRITE( *, 10) PSIZ( I ) 

CONTINUE 
PRINT*, '  ' 

PRINT*, 'THE  PSI  FUNCTIONS  AT  130  ARE:’ 

DO  310  1=1,14 

WRITE(*,10)PSI(I) 

310  CONTINUE 

PRINT*, '  ' 

PRINT*,  'THE  PSI-PRIME  FUNCTIONS  AT  ZERO  ARE:' 
PRINT*,'  ' 

DO  320  1=1,14 

WRITE (*,10)PSIPRZ(I) 

320  CONTINUE 

PRINT*. '  ' 

PRINT*. 'THE  PSI-PRIME  FUNCTIONS  AT  130  ARE:' 
PRINT*,'  ' 

DO  330  1=1,14 

WR ITE ( • , 1 0 ) PS  I  PR ( I ) 

CONTINUE 

ENDIF 


330 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 


STORE  THE  PHIPR(I),  THEPR(I),  PSIPR(I),  PHICRKI),  PHICR2(I), 
THECRl(I),  AND  THECR2(I)  ARRAYS  IN  A  SEQUENTIAL  FILE  CALLED  ’ FRH S ' 
FOR  USE  BY  SHUTBM. 


OPEN ( UN IT= 1 3 , FILE= ' FRH S ' . ACCESS= ' SEQUENTIAL' , STATU S= ' NEW ' ) 


321 

C 


322 


DO  321  1=1  ,14 

WRITE ( 13 , 1 0 ) PH  I  PR ( I ) 
CONTINUE 

DO  322  1=1,14 

WRITE( 13 , 10) THEPR( I) 
CONTINUE 
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c 


£ 


fr 


€1 


4 


4 


« 


DO  323  1-1,14- 

WRITE ( 13 ,10) PS 1(1) 

323  CONTINUE 
C 

DO  324  1-1,14 

WRITE ( 13 ,10) PHICR1 ( I ) 

324  CONTINUE 
C 

DO  325  1-1,14 

WRITE (13 ,10) PH ICR2 ( I ) 

325  CONTINUE 
C 

DO  326  1-1,14 

WRITE (13,10) THECR1 ( I ) 

326  CONTINUE 
C 

DO  327  1-1,14 

WRITE (13,10) THECR2 ( I ) 

327  CONTINUE 
C 

C 

c 

C  STORE  THE  ANGLES  AND  CORRECTION  IN  THE  DATA  FILE  'FPHCRS'  FOR 

C  CLOSE  INSPECTION.  IF  DESIRED: 

C 

c  ***•••*•*•*•*•*•»***«**»•*******•••»•**•**********•*•••**••••»••• 

c 

OPEN  ( UNIT- 9 , FILE- ' FPHCRS ' , ACCESS- ' SEQUENTIAL ' , STATU S= ' NEW ' ) 
WRITE(9. •) 'THE  ALPHAnS  ARE:' 

WRITE ( 9 , * ) '  ' 

DO  311  1-1,14 

WRITE ( 9,10) ALPH A( I ) 

311  CONTINUE 

WRITE ( 9 , • )  '  ' 

WRITE (9,*) 'THE  ALPH A- 40 S  ARE:’ 

WRITE( 9 , • ) '  ' 

DO  312  1-1,14 

WRITE ( 9 , 1 0) ALPHA ( I ) *SN1 

312  CONTINUE 

WRITE ( 9 , • ) ’  ' 

WRITE(9,*) 'THE  PHI  FUNCTIONS  AT  40  ARE:' 

WRITE( 9 , • ) ’  ' 

DO  331  1-1,14 

WRITE ( 9,10)PHICR1(I) 

331  CONTINUE 

WRITE ( 9 , * ) '  ' 

WRITE  ( 9,  •)  'THE  PHI  FUNCTIONS  AT  80  ARE:' 

WRITE( 9  ,  • )  '  ' 

DO  332  1-1,14 

WRITE (9 , 10) PHICR2( I) 


332  CONTINUE 
WRITE ( 9 ,  • )  '  * 

WRITE ( 9. *) ’THE  THETA  FUNCTIONS  AT  40  ARE:' 

WRITE (  9  ,  * ) '  ' 

DO  333  1*1,14 

WRITE (9 ,10) THECR1 (I) 

333  CONTINUE 
WRITE (9,*)'  ' 

WRITE  ( 9,  •)  'THE  THETA  FUNCTIONS  AT  80  ARE:' 

WRITEC9,*)'  ' 

DO  334  1-1,14 

WRITE( 9 , 10 ) THECR2 ( I ) 

334  CONTINUE 
ENDFILE( UNIT-9) 

C 

c 

C  THIS  NEXT  SECTION  CREATES  THE  MATRICES  FR,  6R,  FP,  OP.  6T. 

C  AND  FT.  FR  IS  THE  ROLL  MASS  MATRIX  MODIFYING  TERM  DUE  TO  THE 
C  MOTION  OF  THE  MASSES  ON  EACH  END  OF  THE  BEAM.  FP  AND  6T  ARE 
C  THE  PITCH  AND  TAW  MODIFYING  TERMS  DUE  TO  THE  SAME  THING.  GR. 

C  GP.  AND  FT  (WHICH  ENDS  UP  BEING  ZERO  FOR  THE  FIXED-FREE  MODEL) 
C  ARE  MASS  MATRIX  MODIFYING  TERMS  DUE  TO  THE  COUPLING  BETWEEN 
C  ROLL,  PITCH,  AND  TAW  WHEN  THE  BEAM  IS  DEFORMED. 

C 

c 

DO  340  1-1,14 
DO  350  J-1,14 

FR ( I , J) — ( M1*PH IZ ( J) *PH IZ ( I )  +  M4*PHI( J)*PHI(I) 

1  +  A1*PHIPRZ( J) •PHIPRZ(I)  +  A4  *  PH I PR (J)* PH I PR (I)) 

350  CONTINUE 
340  CONTINUE 
C 

DO  360  1-1,14 
DO  370  J-1.14 

GR( I , J ) -  Z4*(PSI( J) *PHI( I) )  +  D1 *PS IZ ( J ) *PH IPRZ ( I ) 

370  CONTINUE 
360  CONTINUE 
C 

DO  3  80  1-1,14 
DO  390  J-1,14 

FP ( I , J) — ( M1*THETAZ ( J) *THETAZ ( I )  +  M4*THETA( J) *THETA( I ) 

1  +  B1 *THEPRZ( J ) *THEPRZ ( I )  +  B4*THEPR( J) •THEPR( I ) ) 

390  CONTINUE 
3  80  CONTINUE 
C 

DO  400  1-1,14 
DO  410  J-1,14 

GP( I , J) -  T4*(PSI(J) *THETA ( I ) ) 

410  CONTINUE 
400  CONTINUE 


n  n  n  n  n  n  n 


DO  430  J-1,14 

GY(I,J)-( -D1 ) *<PSIZ(I))*( PHIPRZ ( J) ) 

430  CONTINUE 
420  CONTINUE 
C 

DO  440  1-1,14 
DO  450  J-1,14 

FY(I,J)  —  (C1*PSIZ(I)*PSIZ(J)  +  C4*PSI(I)*PSI(J) ) 

450  CONTINUE 
440  CONTINUE 
C 

C  . . . . . . . 

c 

THIS  SECTION  CREATES  THE  ROLL  AND  PITCH  'CORRECTION  MATRICES'. 
THESE  ARE  ACTUALLY  ADDITIONAL  TERMS  ORIGINALLY  IGNORED  IN  THE 
DEVELOPMENT  OF  THE  EQUATIONS  OF  MOTION  FOR  THE  SYSTEARE 
DUE  TO  THE  MOTION  OF  THE  MASSES  IN  THE  PROOF-MASS  ACTUATORS. 


DO  451  1*1,14 
DO  452  J-1,14 

PHICOR( I , J) — M2* (PHICR1 ( I ) *PHICR1 ( J) +  PHICR2 ( I ) *PHICR2 ( J) ) 
THE  COR  (  I,  J)  —M2*  (  THECR1  (I)  *THECR1  (  J)  +  THECR2  (I)  *THECR2(  J)  ) 
452  CONTINUE 
451  CONTINUE 
C 

C  . . . 

c 

OPEN( UNIT- 1 9 . FILE- ' FCORS ' , ACCESS- ' SEQUENTIAL ' , STATUS- ’ NEW ' ) 

16  FORMAT (5(E19.12)) 

17  FORMAT (4(E19.12)) 

WRITE( 19 , • ) '  ' 

WRITE(19,*) 'THE  PHI  CORRECTION  MATRIX  IS:' 

WRITE( 19 , •)  '  ' 

K-l 
L-5 

455  DO  456  1-1,14 

WRITE( 19,16) (PHICOR(I, J) ,J-K,L) 

456  CONTINUE 
WRITE ( 19 , • )  '  ' 

E-K+5 
L-L+5 

IF  (L.NE.15)  GOTO  455 
DO  457  1-1,14 

WRITE ( 19,17) ( PH  I COR ( I , J) , J-11,14) 

457  CONTINUE 
ENDFILE( UNIT-1 9) 


IF  (I.EQ.lll)  THEN 


PRINT*. *  * 

PRINT*. 'THE  U-PHI  ROLL  MATRIX  IS:' 
PRINT*.'  ' 

DO  460  1-1.14 

WRITE ( *.10) (FR(I.J) .J-1,14) 
CONTINUE 
PRINT*. '  ' 

PRINT*. 'THE  U-PSI  ROLL  MATRIX  IS:' 
PRINT*.'  ' 

DO  470  1-1.14 

WRITE (*,10) (GR(I.J) .J-1,14) 
CONTINUE 


PRINT*. '  ' 

PRINT*.  'THE  U- THETA  PITCH  MATRIX  IS 
PRINT*,'  ' 

DO  480  1-1.14 

WRITE(*.10) (FP(I.J) .J-1,14) 
CONTINUE 
PRINT*. '  ' 

PRINT*. ' THE  U-PSI  PITCH  MATRIX  IS:' 
PRINT*.'  ' 

DO  490  1-1,14 

WRITE ( *,10) (GP(I,J) .J-1,14) 
CONTINUE 

PRINT*. 'THE  U-PHI  TAW  MATRIX  IS:’ 
PRINT*,'  ' 

DO  500  1-1,14 

WRITE(*,10) (GY(I. J) , J-l .14) 
CONTINUE 
PRINT*. '  ' 


PRINT*. 'THE  U-PSI  TAW  MATRIX  IS:' 
PRINT*,'  ' 

DO  510  1-1,14 

WRITE(*. 10) (FT(I.J) .J-1,14) 
CONTINUE 


ENDIF 


THIS  SECTION  CREATES  THE  MASS  MATRIX  'M'  AND  OMEGA- SQUARED 
MATRIX  'O'.  WHICH  WILL  HAVE  AN  EIGENVALUE  PROBLEM  DONE  ON  THEM. 


DO  520  I-l.V 
DO  530  J-l.V 

IF  (I.EQ.J)  THEN 

M(I. J)-1-FR(I. J)- PHI COR (I, J) 

ELSE 

M(I,J)— FR(I.J) -PH I COR ( I , J) 
ENDIF 
CONTINUE 
CONTINUE 

DO  540  I-l.V 
DO  550  J-W.X 
M( I , J) *0 
CONTINUE 
CONTINUE 

DO  560  I-l.V 
DO  570  J-Y.Z 
L-J-X 

M(I,  J)—  GR(I.L) 

CONTINUE 

CONTINUE 


DO  5  80  I-W.X 
X-  I-V 

DO  590  J-l.V 
M ( I , J) -0 
CONTINUE 
CONTINUE 

DO  600  I-W.X 
X-  I-V 

DO  610  J-W.X 
L-J-V 

IF  (I.EQ.J)  THEN 
M(I,J)-1-FP(K.L)  -  THE  COR  (  X.  L) 
ELSE 

M( I , J) — FP (X, L) -THEC0R( X, L) 
ENDIF 
CONTINUE 
CONTINUE 

DO  620  I-W.X 
X-  I-V 

DO  630  J-T.Z 


MCI.  J)—  GP(I.L) 

630  CONTINUE 
620  CONTINUE 
C 
C 

DO  640  I-Y.Z 
I-I-X 

DO  650  J-l.V 

M(I, J)«-6Y(I, J) 

650  CONTINUE 
640  CONTINUE 
C 

DO  660  I-Y.Z 
DO  670  J-W.X 
«(I. J)-0 
670  CONTINUE 
660  CONTINUE 
C 

DO  6  80  I-Y.Z 
I-I-X 

DO  690  J- Y, Z 
L-J-X 

IF  (I.EQ.J)  THEN 
M(I.  J)»1-FY(I,L) 

ELSE 

H(I,  J)—  FY(K.L) 

ENDIF 

690  CONTINUE 
680  CONTINUE 
C 

C  WRITE  THE  M  MATRIX 
C 

IF  (I.EQ.lll)  THEN 
PRINT*. 'THE  M  MATRIX  IS:' 

PRINT*. '  ' 

DO  700  I-l.Z 
DO  710  J-l.Z 

WRITE ( *,10)M(I.J) 

710  CONTINUE 
700  CONTINUE 

ENDIF 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  PUT  THE  M  MATRIX  IN  A  FILE  CALLED  FMASTM 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

OPEN( UN  IT-8,  FILE-' FMASTM' , ACCESS- ' SEQUENTIAL' . STATU S- ' NEW ' ) 
DO  711  I-l.Z 
DO  712  J-l.Z 

WRITE ( 8,10)M(I,J) 

712  CONTINUE 

711  CONTINUE 
ENDFILE( UNIT-8) 


*• 


ccccccccccccccccccccccccccccccccc 

C  CREATE  THE  OMEGA- SQUARED  MATRIX 

ccccccccccccccccccccccccccccccccc 

c 

DO  720  I-l.V 
DO  730  J-l.Z 

IF  (I.EQ.J)  THEN 
0(  I.  J)-OMESQR(I) 

ELSE 

0(1. J)-0 

ENDIF 

730  CONTINUE 
720  CONTINUE 
C 

DO  740  I-W.X 
I-  I-V 

DO  750  J-l.Z 

IF  (I.EQ.J)  THEN 
0(1. J)-OMESQP(K) 

ELSE 

0(1. J)-0 
ENDIF 
750  CONTINUE 
740  CONTINUE 
C 

DO  760  I-T.Z 
I-I-X 

DO  770  J-l.Z 

IF  (I.EQ.J)  THEN 
0(1. J)-OMESQT(K) 

ELSE 

0(1. J)-0 

ENDIF 

770  CONTINUE 

760  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  WRITE  THE  0  MATRIX  INTO  A  FILE  CALLED  FMASTX 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

OPEN< UNIT-11 .FILE- 'FMASTX' , ACCESS- ' SEQUENTIAL • , STATU S- *  NEW  * ) 
C 

DO  761  I-l.Z 
DO  762  J-l.Z 

WRITEdl,  10)0(1,  J) 

762  CONTINUE 

761  CONTINUE 
ENDFILE( UNIT-11 ) 


DO  7  80  I-l.Z 
DO  790  J-l.Z 


C 


AA  (  I ,  J )  =  REAL  ( 0  ( I ,  J  )  ) 
BB ( I , J ) = REAL (  M  ( I ,  J )  ) 
CONTINUE 
CONTINUE 


THIS  SECTION  BUILDS  THE  [FB]  MATRIX  ' FBMAT' ,  WHICH  IS  WHAT  IS 
LEFT  OVER  ON  THE  RIGHT  HAND  SIDE  (OTHER  THAN  THE  INVERSE  TRANS¬ 
FORMATION  MATRIX)  AFTER  THE  CONTROL  VECTOR  [FC]  IS  BROKEN  OUT. 
THE  RIGHT  HAND  SIDE  LOOKS  LIKE:  RHS  =  [E][FB]FC. 


1010 

1000 

C 


DO  1000  1=1 ,45 
DO  1010  J=l,10 
FB  MAT ( I , J ) =0 
CONTINUE 
CONTINUE 

FB  MAT (1,2)= SN1 
FBMAT ( 1 , 4)  =  SN2 
FBMAT( 1.5) =1000.0 
FBMAT(1.6)=1000.0 
FBMAT (2,1)=- SN1 
FBMAT  (2,3)=-  SN2 
FBMAT (2,7)=1000.0 
FB  MAT  (  2 , 8)  =1000 .0 
FBMAT( 3,9)=1000.0 
FBMAT (3 ,10) *1000 .0 


DO  1020  1=4.17 
K-I-3 

FBMAT (1,1) =-PHICRl ( K) 

FB  MAT ( 1 , 3 )  “-  PH  I CR2  ( K ) 

FBMAT (I, 6) =-PHIPR(K)*1000 .0 
1020  CONTINUE 


DO  1030  1=18,31 
K=  1-17 

FBMAT( I . 2 ) — THECR1 ( K) 

FBMAT( 1,4) — THECR2(K) 

FBMAT (1,8) *-THEPR (K)*1000.0 
CONTINUE 


1030 


32,35 


DO  1040  1=32, 
K= 1-31 
FBMAT (1,10) 


!PS I ( K ) *  1000 . 0 


1040  CONTINUE 
C 

DO  1050  1*1,45 
DO  1060  J-1,10 

FBMATK  I,  J)  =  REAL (FB MAT ( I,  J)  ) 

1060  CONTINUE 
1050  CONTINUE 
C 

c 

C  THIS  SECTION  CREATES  THE  P-MATRIX,  WHICH  IS  FROM  THE  EQUATION 

C  T*P*X.  THE  P-MATRIX  WILL  BE  MULTIPLIED  BY  THE  MATRIX  OF  EIGEN- 
C  VALUES  IN  ' SHUTBM'  TO  PUT  THE  STATE  VECTOR  IN  MODAL  COORDINATES. 
C 

c 

c 

DO  1070  1*1,7 
DO  1080  J*1 , 45 
PMAT( I, J) *0 
10  80  CONTINUE 
1070  CONTINUE 
C 

PMAT ( 1 , 1 ) *1 . 0 
PMAT( 2 , 2 ) *1 . 0 
PMAT ( 3 , 3 ) *1  .0 
C 

DO  1090  J-4.17 
I-J-3 

PMAT ( 4 , J) =PHICR1 (X) 

PMAT (6,J)=PHI CR2 ( K ) 

1090  CONTINUE 
C 

DO  1100  J*1 8 , 3 1 
X* J-17 

PMAT( 5 , J) *THECR1 (I) 

PMAT ( 7 , J ) *  THECR2 ( K ) 

1100  CONTINUE 
C 

DO  1110  1*1,7 
DO  1111  J=1 , 45 

PMATX(I, J)*REAL(PMAT(I, J) ) 

1111  CONTINUE 
1110  CONTINUE 
C 

c 

C  STORE  THE  FBMATX  AND  THE  PMATX  IN  A  FILE  CALLED  * FBMATX ' : 

C  M*MM****MMMMMMM*M*MM*******MMMM**M«tMM«*M*M 

c 

c 

OPEN{ UNIT-12,  FILE* 'FBMATX*  , ACCESS* ’ SEQUENTIAL  * , STATU S- *  NEW • ) 

C 
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2  FORMAT ( E14 . 7) 

DO  1120  1*1.45 

WRITE (12,12)  (FBMATXCI, J)  .J-1,10) 
120  CONTINUE 

DO  1130  1=1.7 

WRITE  (12. 12)  (  PMATX(  I ,  J)  .  J=l,45) 
130  CONTINUE 


ENDFILECUNIT-12) 


C  THIS  SECTION  SOLVES  THE  EIGENVALUE  PROBLEM  A*X  *  LAMB DA *B  *X, 

C  WHERE  A  IS  THE  'O'  MATRIX  AND  B  IS  THE  'M'  MATRIX.  IT  USES  THE 
C  IMSL  ROUTINE  EIGZF.  WHICH  RETURNS  THE  EIGENVALUES  ( EIGENV)  AND 
C  EIGENVECTORS  (VEC)  . 


IJ  OB*  2 

CALL  EIGZF  (AA. IA. BB, IB, N. IJOB, ALFA, BBETA. VEC. IZ, WK, IER) 


DO  800  1=1. Q 

IF  (BBETA(I) .NE.O)  THEN 

EIG ENV ( I )  =  ALFA ( I ) / BBETA ( I ) 
ELSE 

EIGENV (I)=999999.9 
ENDIF 

800  CONTINUE 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  SORT  THE  EIGENVALUES  AND  EIGENVECTORS 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 

IF  (I.NE.lll)  GOTO  861 

U=Q-1 

810  IF  (U.NE.O)  THEN 
DO  820  1=1, U 
K=I+1 

COMF1 =REAL( EIGENV (I) ) 


C0MP2 - REAL (E 16 ENV ( I ) ) 

IF  (C0MP1.LT. C0NP2)  THEN 
TEMP(I)=EIGENV(I) 

m  E 16 ENV  ( I )  -EIGENV  ( K) 

EI6  ENV ( K) » TEMP ( I ) 

DO  830  J-l.Q 

VECT( J, I)-VEC( J, I) 

VEC( J, I) =VEC( J,  K) 

VEC( J,K)=VECT( J. I) 
m  830  CONTINUE 

ENDIF 

820  CONTINUE 

u=  u-i 

GOTO  810 
ENDIF 

•  c 

861  CONTINUE 
C 
C 
C 

^  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  MAKE  THE  DOUBLE  PRECISION  NUMBERS  INTO  REAL 
C  CALCULATE  THE  FREQUENCIES  IN  RADIANS  PER  SECOND  (OMEGA) 

C  CALCULATE  THE  FREQUENCIES  IN  CTCLES  PER  SECOND  (HERTZ) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

•  DO  840  1*1,0 

DO  850  J-l.Q 

VREAL ( I , J) =REAL ( VEC ( I , J) ) 

850  CONTINUE 
840  CONTINUE 
C 

•  DO  860  1=1, Q 

OMEG  A  ( I )  =  SQ  RT  ( E IG  ENV  ( I )  ) 

HERTZ ( I ) =OMEG A ( I ) / ( 2  *ACOS( -1 . 0 ) ) 

860  CONTINUE 

C 

C 

fc  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C  STORE  THE  INFORMATION  IN  FILES  ' FFREQQ'  AND  'FIGVCT'. 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c 

OPEN( UNIT-15 , FILE- *  FFREQQ 1 , ACCESS- ’ SEQUENTIAL' , STATUS- ’ NEW • ) 

•  7  FORMAT (A2 4) 

WRITE( 15,7) * 

WRITEU5 , 7)  ’  THE  EIGENVALUES  ARE: 

WRITE( 15,7) * 

DO  870  1-1,42 

WRITE ( 15,12) EIGENV ( I ) 

•  870  CONTINUE 
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non 


ENDFILE( 0NIT=1 5 ) 


C 

OPEN ( UNIT= 10 , FILE= ' FIGVCT' . ACCESS= ' SEQUENTIAL ' , STATU S= ' NEW ' ) 
C 

11  FORMAT ( 6 ( El 4 . 7 ) ) 

C 

WRITE ( 10, •) 'THE  UNSORTED  EIGENVECTOR  MATRIX  IS:' 

WRITE(10,*)'  ' 

K=1 

L=6 

950  DO  940  1=1,42 

WRITE (10, 11) (VREAL(I, J) , J=K,L) 

940  CONTINUE 

WRITE (10 , ♦ )  '  ' 

K=K+6 

L=L+6 

IF  (L.NE.48)  GOTO  950 
ENDFILE( UNIT=10) 


C 


PROGRAM  SHUTBM 


C 

C  THIS  PROGRAM  SOLVES  AN  EIGENVALUE  PROBLEM  OF  THE  FORM 

C  A*X  =  LAMB  DA*B  *X .  THE  A  MATRIX  IS  A  45X45  MATRIX  MADE  UP 
C  OF  THE  42X42  MATRIX  'O'  FROM  FMASTK.  THE  UPPER  LEFT  CORNER 
C  HAS  THREE  ZEROS  ON  THE  DIAGONAL.  THE  B  MATRIX  IS  A  45X45 
C  MADE  OF  FOUR  PARTS.  THE  UPPER  LEFT  3X3  IS  THE  MOMENT  OF  INERTIA 
C  MATRIX  FROM  THE  NASA  PAPER.  THE  LOVER  RIGHT  42X42  IS  THE  M 
C  MATRIX  FROM  FMASTM.  THE  UPPER  RIGHT  3X42  (THE  LOWER  LEFT  42X3 
C  IS  THE  TRANSPOSE)  IS  A  CORRECTION  MATRIX  WHICH  TAKES  INTO  AC- 
C  COUNT  THE  ROTATION  OF  THE  SHUTTLE- BEAM- REFLECTOR  SYSTEM.  IT 
C  IS  FILLED  WITH  ZEROS  EXCEPT  FOR  THE  LOWER  LEFT  CORNER.  WHICH 
C  CONTAINS  THE  SZ'S. 

C 

DOUBLE  PRECISION  DALFAL ( 14 ) . DPM ( 42 . 42 ) , DPO ( 42 , 42 ) 

DOUBLE  PRECISION  DPH IPR ( 1 4) . DTHEPR ( 1 4 ) , DPSI ( 1 4 ) , DHI CR1 ( 1 4 ) 
DOUBLE  PRECISION  DHICR2 ( 14 )  . DHECR1 ( 14 ) , DHECR2 ( 14 ) 

REAL  BETA( 45 ) , ALPHAL( 14) . M ( 42 , 42 ) , 0 ( 42 . 42 ) , A( 45 , 4 5 ) , B( 45 . 45 ) 

REAL  WK( 2888 ) ,SZ(14) ,SZMAT(3,42) ,INRTIA(3,3) ,EIGENV(45) 

REAL  EIGVEC(45 ,45)  .  MTILDA(  45 , 45  )  »  FBMATX(45 . 10)  ,  PRODK45 .45) 

REAL  INVRSE( 45,45) , PROD2 (45,45) , PROD3 (45,45) , KTILDA ( 45 , 4  5 ) 

REAL  DTILDA( 45 , 45 ) , PHIPR( 14) , THEPR( 14) , PS 1(14)  ,PHICR1(14) 

REAL  PH ICR2 ( 14 ) , THECR1 ( 14 ) ,THECR2(14) ,DMAT(45,10) , A4 , B4 . C4 
REAL  CMAT(7 .45) , PMAT(7 .45) 

C 

COMPLEX  ALFA ( 45 ) ,DIGENV(45) , DIGVEC ( 45 , 4 5 ) ,OMEGA(45) 

INTEGER  I.J.K.L.IDGT  ,  LL,  MM,  NN,  IC 
INTEGER  IA.  IB,  N.  IJOB,  IZ,  IER 
C 

PARAMETER  ( A4-496 9 . 0 . B4=4 96 9 .0 , C4-993 8 . 0 ) 

C 

c 

C  READ  IN  THE  INFORMATION  FROM  ' FDATA2 '  (THE  ANGLES), 

C  'FMASTM'  (THE  MASS  MATRIX).  AND  'FMASTK'  (THE  STIFFNESS  MATRIX) 

C 

C  . . . . . 

C 

OPEN ( UNIT* 7 . FILE= ' FDATA2 ' , ACCESS* ' SEQUENTIAL' , STATUS* ’ OLD' ) 
OPEN(UNIT=8, FILE*’ FMASTM' , ACCESS* ' SEQUENTIAL' , STATU S* ' OLD ' ) 
OPEN(UNIT=9. FILE* 'FMASTK' .ACCESS*’ SEQUENTIAL' , STATUS* 'OLD' ) 
OPEN( UNIT=12 , FILE* ' FRHS ' , ACCESS* ' SEQUENTIAL’ , STATU S* ' OLD ' ) 

OPEN ( UNIT=11 .FILE* ' FBMATX ' , ACCESS* ' SEQUENTIAL' , STATUS* ' OLD' ) 
REWIND  (UN  IT*  7) 

REWIND  (UNIT*  8) 

REWIND  (UN  IT*  9) 

REWIND ( UNIT-12) 

REW  IND (  UNIT*  11) 

C 

10  FORMAT(E31  .24) 

11  FORMAT (El  4 .7) 

13  FORMAT (10(E11 .4) ) 


FORMAT (  4  (  3  (  4X, E14 .7) . /) , /) 
FORMAT (S(E14.7) ) 

DO  20  1*1.14 

READ (7,10) DAL  FAL ( I ) 

ALPBAL ( I ) *  REAL ( DALF AL ( I ) ) 
CONTINUE 

DO  30  1-1,42 
DO  40  J-1.42 

READ (8,10) DPM ( I , J ) 

M ( I , J ) *  REAL ( DPM ( I , J ) ) 
READ (9,10) DPO ( I , J ) 

0(1, J)=REAL(DPO(I, J) ) 
CONTINUE 
CONTINUE 

DO  41  1-1,14 

READ (12,10) DPHIPR ( I ) 

PHIPR ( I ) *REAL( DPHIPR ( I ) ) 
CONTINUE 

DO  42  1*1,14 

READ (12.10) DTHEPR ( I ) 

THEPR ( I ) - REAL ( DTHEPR ( I ) ) 
CONTINUE 

DO  43  1*1,14 

READ ( 12 . 10 ) DPS I (I) 

PSI ( I ) -REAL( DPS  I ( I ) ) 
CONTINUE 

DO  44  1-1,14 

READ (12,10) DH ICR1 ( I ) 
PHICR1 ( I ) -REAL ( DH ICR1 ( I ) ) 
CONTINUE 

DO  43  1-1,14 

READ ( 1 2 , 1 0 ) DH I CR2 ( I ) 

PH  I CR2 ( I ) *REAL( DH I CR2 ( I ) ) 
CONTINUE 

DO  46  1-1,14 

READ (12,10) DHECR1 ( I ) 
THECR1 ( I ) -REAL( DHECR1 ( I ) ) 
CONTINUE 

DO  47  1-1,14 

READ( 12,10) DHECR2 ( I ) 
THECR2 ( I ) -REAL( DHECR2 ( I ) ) 
CONTINUE 


DO  48  1-1,45 

READ (11,11) (FBMATX(I, J) , J-1,10) 

48  CONTINUE 
C 

DO  49  1-1,7 

READ (11,11) (PMAT(I, J) , J-1.45) 

49  CONTINUE 
C 

ENDFILE( UNIT-7) 

ENDFILE( UNIT-8) 

ENDFILE( UNIT-9) 

ENDFILEC UNIT-12) 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCC 
C  FORM  THE  INRTIA  MATRIX: 

CCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

INRTIA (1,1)- 97 5423  .0 
INRTIA( 1 , 2 ) -0 
INRTIAd, 3)— 145393  .0 
INRTIA ( 2, l)-0 
INRTIA (2. 2) -6  85  90  80  .0 
INRTIA( 2 , 3 ) -0 
INRTIA (3.1)— 1453  93  .0 
INRTIAd  ,2)  =0 
INRTIA (3, 3) -7086601  .0 
C 

CCCCCCCCCCCCCCCCCCCCCC 
C  CALCULATE  THE  SZ'S: 

CCCCCCCCCCCCCCCCCCCCCC 

C 

DO  50  1-1.14 

SZ(I)-(2*130*SQRT( .09556*130) ) / ( ALPHAL(I) **2) 

50  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCC 
C  BUILD  THE  SZMAT  MATRIX: 

CCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

DO  60  1-1,3 
DO  70  J-1,42 
SZMAT ( I . J) -0 
70  CONTINUE 

60  CONTINUE 

C 
C 
C 
C 

DO  80  J-1,14 

SZMAT(1. J)-(A4*PHIPR( J) ) 

SZMAT(  2  ,  J)  -+  SZ(  J) 

80  CONTINUE 


DO  81  J-15,28 
K-J-14 

SZMAT (  2  ,  J ) -B4 *THEPS ( K ) 
SZ  MAT  ( 1 ,  J )  —  SZ  ( K ) 
CONTINUE 

DO  82  J-29,42 
I-J-28 

SZMAT(3.J)-C4*PSI(K) 

CONTINUE 


THIS  SECTION  FORMS  THE  LARGE  'A'  MATRIX,  WHICH  CONSISTS  OF 
THE  ONEGA- SQUARES  ON  THE  DIAGONAL  WITH  ZEROS  AS  THE  FIRST  THREE 
OMEGAS.  IT  THEN  BUILDS  THE  '  B'  MATRIX.  WHICH  IS  MADE  UP  OF 
THE  INERTIA  MATRIX.  THE  MASS  MATRIX.  AND  THE  SZ-MATRIX  (SZMAT) 
IT  WILL  LOOK  LIKE: 
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DO  100  I-1.4S 
DO  110  J-1,45 
A(I. J)-0 
CONTINUE 
CONTINUE 


DO  120  1-4. 45 


I 


DO  130  J-4.45 
L-J-3 

A(I,  J)-O(K.L) 

130  CONTINUE 
120  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  CREATE  THE  LARGE  B  MATRIX: 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 
C 

DO  140  1-1.3 
DO  130  J-1.3 

B ( I , J ) - INRTIA ( I . J ) 

150  CONTINUE 
140  CONTINUE 
C 

DO  160  1-1.3 
DO  170  J-4.45 
L-J-3 

B(I.  J)-SZMAT(I.L) 

170  CONTINUE 
160  CONTINUE 
C 

DO  180  1-4.45 
X-  1-3 

DO  190  J-1.3 

B(I,J)-SZ  MAT  (  J .  K ) 

190  CONTINUE 
180  CONTINUE 
C 

DO  200  1-4.45 
I- 1-3 

DO  210  J-4.45 
L-J-3 

B(I, J)-M(I.L) 

210  CONTINUE 
200  CONTINUE 
C 

cccccccccccccccccccccccccccccccccccccc 

C  CHECK  THE  INRTIA  AND  SZMAT  MATRICES: 

cccccccccccccccccccccccccccccccccccccc 

c 

IF  (I.EQ.lll)  THEN 
PRINT*. 'THE  INERTIA  MATRIX  IS:' 
PRINT*. '  ' 

DO  220  1-1.3 

VRITE( *.11) ( INRTIA ( I ,  J ) .J-1.3) 
220  CONTINUE 


u  u  u  u 


19  FORMAT (3(E14.7) ) 

C 

PRINT*. *  ' 

PRINT*. ' THE  SZ  TRANSPOSE  MATRIX  IS:' 

PRINT*. '  * 

DO  221  1-1,42 

WRITE(  *,19)  (  SZ  MAT  ( J ,  D.J-1.3) 

221  CONTINUE 
C 

ENDIF 

ccccccccccccccccccccccccccccccccccccccccccccccccccccc 


THIS  SECTION  SOLVES  THE  EIGENVALUE  PROBLEM  A*X  -  LAMBDA*B*X, 
GIVING  BACK  THE  EIGENVALUES  AND  EIGENVECTORS  OF  THE  SYSTEM. 


CALL  SORT(EIGENV,  EIGVEC,  45 ) 

C 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  CHECK  THE  EIGENVALUE  AND  EIGENVECTORS : 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

OPEN( UNIT-10,  FILE- 'FIG ENV ' , ACCESS* ' SEQUENTIAL ' . STATU S« ' NEW ' ) 
C 

17  FORMAT (5(E14.7)) 

DO  260  1*1 ,45 

WRITE ( 10,17) EIGENV ( I ) 

260  CONTINUE 

WRITE (10,*)'  ' 

C 

K*1 

L*5 

261  DO  262  1*1 ,45 

WRITE (10, 17) (EIGVEC ( I , J) . J-K.L) 

262  CONTINUE 

WRITE( 10  ,  • )  '  ' 

1*1+5 

L»L+5 

IF  (L.NE.50)  GOTO  261 
ENDFILE( UNIT*10) 

C 

C 

C 

998  LL*45 

MM*  4  5 
NN*45 
IC-45 
IDGT-5 
C 

. . . 

c 

c 

c 

c 

c 

c 

ccccccccccccccccccccccccccccccccccccccccccc 

C  MULTIPLY  B •TRANSFORMATION  MATRIX: 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL  VMULFF( B, EIGVEC. LL. MM, NN, IA, IB, PROD1, IC, IER) 

C 

PRINT*, ' THE  IER  FOR  VMULFF  (B, EIGVEC)  IS:  '.IER 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  INVERT  THIS  PRODUCT 


ccccccccccccccccccccccccccccccccccccccccc 

c 

CALL  LINV2F  ( PE0D1 ,  N,  I  A,  INVRSE,  IDGT,  WI,  IER) 

PRINT* , *  THE  IER  FOR  LINV2F  (PROD1)  IS:  '.IER 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  MULTIPLY  THE  INVERSE  AND  B: 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 

CALL  VMULFF  ( INVRSE. B. LL. MM. NN. IA. IB. PROD2 . IC, IER) 

C 

PRINT*. 'THE  IER  FOR  VMULFF  ( INVERSE, B)  IS:  '.IER 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  MULTIPLY  THIS  RESULT  BY  THE  TRANSFORMATION  MATRIX: 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

CALL  VMULFF(PROD2 , EIGVEC. LL. MM, NN. IA, IB.MTILDA. IC. IER) 

C 

PRINT*, 'THE  IER  FOR  VMULFF  ( PROD2 , EIGVEC)  IS:  '.IER 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  PREMULTIPLY  THE  INVERSE  TRANSFORMATION  MATRIX  BY  A: 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

CALL  VMULFF( INVRSE, A. LL, MM, NN, IA. IB.PROD3.IC, IER) 

C 

PRINT*. ' THE  IER  FOR  VMULFF  (INVERSE. A)  IS:  '.IER 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  MULTIPLY  THIS  RESULT  BY  THE  TRANSFORMATION  MATRIX 

cccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

CALL  VMULFF(PROD3, EIGVEC. LL. MM, NN, IA. IB.KTILDA, IC. IER) 

C 

PRINT*. 'THE  IER  FOR  VMULFF  (PROD3 , EIGVEC)  IS:  '.IER 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  MULTIPLY  INVRSE  AND  FBMATX  TO  GET  DMAT,  WHICH  WILL  BE  THE  RHS 
C  MATRIX  IN  ' FCNTRL ' . 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

CALL  VMULFF( INVRSE. FBMATX, LL, MM, 10. IA, IB. DMAT, IC, IER) 

C 

PRINT*, 'THE  IER  FOR  VMULFF  ( INVRSE, FBMATX)  IS:  '.IER 
C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  MULTIPLY  PMAT  AND  EIGVEC  TO  GET  CMAT,  WHICH  WILL  BE  READ  INTO 
C  ACOSS  DIRECTLY: 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

2222  CALL  VMULFF  (PMAT, EIGVEC. 7 . 43 . 45 ,  7, 45. CMAT.  7, IER) 


c 

PRINT*, 'THE  IER  FOR  VMULFF  ( PMAT, EIGVEC )  IS:  ' . IER 
C 
C 

c 

C  PUT  MTILDA  AND  KTILDA  THROUGH  A  FILTER  TO  TARE  OUT  THE  VERT 
C  SMALL  NUMBERS  (SET  THEM  TO  ZERO) 

C 

c 

DO  510  1=1,45 
DO  520  J=1 , 45 

IF  (MTILDA(I.  J)  .LT.1E-4)  THEN 
MTILDAd,  J)=0 
ENDIF 

IF  (KTILDAd,  J)  .  LT.1E-4)  THEN 
KTILDAd,  J)=0 
ENDIF 

520  CONTINUE 

510  CONTINUE 

C 

C 

C 

C 

c  ***••*••••»•••••••••••••***••*•**•»•••»•••••*»»•••••••••••»••»»*» 

c 

C  THIS  SECTION  CREATES  THE  DAMPING  MATRIX  ' DTILDA* .  IT  TAKES 
C  THE  SQUARES  OF  THE  EIGENVALUES.  WHICH  ARE  ON  THE  DIAGONAL  OF 
C  'KTILDA'.  AND  FORMS  2  *ZETA*OMEG A,  WHERE  ZETA  IS  A  CONSTANT  .003, 
C  AND  OMEGA  IS  THE  SQUARE  ROOT  OF  EACH  EIGENVALUE.  THE  RESULT 
C  IS  A  DIAGONAL  MATRIX  CORRESPONDING  TO  THE  DIAGONAL  STIFFNESS 
C  MATRIX  (KTILDA). 

C 

C  . . . . . . 

C 

c 

DO  530  1=1,45 
DO  540  J=1 ,45 

IF  (I.EQ.J)  THEN 

DTILDA (I, J) =2.0 *.003 *SQRT( KTILDA (I, J)) 

PT  QP 

DTILDA ( I , J) =0 
ENDIF 

540  CONTINUE 
530  CONTINUE 
C 

OPEN ( UNIT- 15. FILE- ' FMTILD' , ACCESS- ' SEQUENTIAL' . STATUS-' NEW ' ) 
OPEN( UNIT-17 . FILE- ' FKTILD* . ACCESS- 'SEQUENTIAL' , STATUS- ' NEW ' ) 
OPEN ( UNIT-1 8 , FILE- ' FDTILD' , ACCESS- ' SEQUENTIAL ' , STATUS- ' NEW ' ) 
OPEN( UNIT- 16 , FILE- ' FTRANS ' . ACCESS- 'SEQUENTIAL' , STATUS- ' NEW ' ) 
OPEN ( UN IT- 13 . FILE- ' FDMATX ' . ACCESS- ' SEQUENTIAL' , STATUS- ' NEW ' ) 
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4 


4 
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C 

K=1 
L=  5 

WRITE ( 15#*) ’ THE  DIAGONAL  MASS  MATRIX-UNSORTED- IS : ' 
WRITE(17»*) 'THE  DIAGONAL  STIFFNESS  MATRIX-UNSORTED- IS : ' 
WRITE(18, •) 'THE  DIAGONAL  DAMPING  MATRIX- DNSORTED- IS : ' 
WRITE ( 15  ,  • )  '  ' 

WRITE ( 17  »  * ) '  ' 

WRITE( 18,  *)  '  ' 

499  DO  500  1=1,45 

WRITE (15,16) (MTILDA(I. J) ,J=K,L) 

WRITE( 17,16) ( KTILDA ( I, J) , J=I, L) 

WRITE (18,16) (DTILDA(I, J) , J-K.L) 

500  CONTINUE 
WRITE (15 , * ) '  ' 

WRITE( 17 , •) '  ' 

WRITE (18.*)'  ' 

1=1+5 

L=L+5 

IF  (L.NE.50)  GOTO  499 
C 

DO  600  1=1,45 

WRITE  ( 16 , 1 6 )  (INVRSEd,  J)  ,J=1,45) 

600  CONTINUE 
C 

WRITE (13 , • )  1,  1.  1 
WRITE ( 13 , •)  1.0 

WRITE ( 13 , • )  6.  4,  2  .10.  7  ,0.003 
DO  601  1=1,12 

WRITE(  13,14)  (EIGVECd.  J)  ,  J-1,12) 

601  CONTINUE 

WRITE ( 13 , '(//)  ') 

C 

C 

DO  610  1-1,12 

WRITE(13 .13) (DMAT(I, J) . J=l,10) 

610  CONTINUE 

WRITE( 13 ,'(///)') 

C 

C  WRITE  IN  THE  FIRST  12  ROWS  OF  CMAT  TRANSPOSED,  WHICH  IS  THE 
C  PROPER  FORM  FOR  ACOSS  TO  HANDLE 
C 

DO  611  1-1,12 

WRITE(13 .13)  (CMATd.  I)  .  J-l  ,7) 

611  CONTINUE 

WRITE ( 13 , '(///)  ') 

C 

C 

2218  FORMAT( 2F3 . 1 ) 

2219  FORMAT( 12F3 . 1) 

DO  620  1-1,12 


si 


4 

r 
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f  RITE(  13 . 11 )  SQRT(  KTILDA  (  I,  I)  ) 

620  CONTINUE 
C 

VRITE( 13  »*<//)  ') 

DO  2200  1=1.4 

WRITE (13 ,2219) (0.0. J-1.12) 

2200  CONTINUE 
WRITE(13, ' (//) ') 

DO  2201  1=1,12 

WRITE (13 ,2218) (0.0. J- 1.2) 

2201  CONTINUE 
C 

ENDFILE( UNIT=1 5 ) 

ENDFILE( UNIT-17) 

ENDFILE( UNIT-18) 

ENDFILE( UNIT-16) 

ENDFILE( UNIT-13 ) 

C 

999  END 
C 

SUBROUTINE  SORT(A.B.D) 

C 

REAL  A( 45 ) ,B(45.45) ,TEMP(45) ,VECT(45,45) 
INTEGER  I.J.K.D.U 
C 

U-D-l 

20  IF  (U.NE.O)  THEN 
DO  30  1=1. U 
1=1+1 

IF  (A(I) .LT.A(I))  THEN 
TEMP ( I ) =A ( I ) 

A( I ) =A(K) 

A ( K) -TEMP ( I ) 

DO  40  J-1.45 

V  E  CT ( J , I)=B(J, I) 

B(  J. I) =B( J,  K) 

B( J.I)=VECT( J. I) 

40  CONTINUE 

ENDIF 

30  CONTINUE 

U-U-l 
GO  TO  20 
ENDIF 
END 
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